Astronomy & Astrophysics manuscript no. ccs 


February 1, 2008 


(DOI: will be inserted by hand later) 





Representations of celestial coordinates in FITS 

Mark R. Calabretta (1) and Eric W. Greisen (2) 

1 Australia Telescope National Facility, P.O. Box 76, Epping, NSW 1710, Australia 

2 National Radio Astronomy Observatory, P.O. Box O, Socorro, NM, USA 87801-0387 



O 
O 

(N 



OS 



Draft dated 17 Jul 2002 

Abstract. In Paper I, Greisen & Calabretta (2002) describe a generalized method for assigning physical coordinates to FITS 
image pixels. This paper implements this method for all spherical map projections likely to be of interest in astronomy. The 
new methods encompass existing informal FITS spherical coordinate conventions and translations from them are described. 
Detailed examples of header interpretation and construction are given. 
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1. Introduction 

This paper is the second in a series which establishes conven- 
tions by which world coordinates may be associated with FITS 
(Hanisch et al., 2001) image, random groups, and table data. 
Paper I (Greisen & Calabretta, 2002) lays the groundwork by 
developing general constructs and related FITS header key- 
words and the rules for their usage in recording coordinate in- 
formation. In Paper III, Greisen et al. (2002) apply these meth- 
ods to spectral coordinates. Paper IV (Calabretta et al., 2002) 
extends the formalism to deal with general distortions of the co- 
ordinate grid. This paper, Paper II, addresses the specific prob- 
lem of describing celestial coordinates in a two-dimensional 
projection of the sky. As such it generalizes the informal but 
widely used conventions established by Greisen (1983, 1986) 
for the Astronomical Image Processing System, hereinafter re- 
ferred to as the AIPS convention. 

Paper I describes the computation of world coordinates as 
a multi-step process. Pixel coordinates are linearly transformed 
to intermediate world coordinates that in the final step are trans- 
formed into the required world coordinates. 

In this paper we associate particular elements of the inter- 
mediate world coordinates with Cartesian coordinates in the 
plane of the spherical projection. Figure 1, by analogy with 
Fig. 1 of Paper I, focuses on the transformation as it applies 
to these projection plane coordinates. The final step is here di- 
vided into two sub-steps, a spherical projection defined in terms 
of a convenient coordinate system which we refer to as native 
spherical coordinates, followed by a spherical rotation of these 
native coordinates to the required celestial coordinate system. 
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Fig. 1. Conversion of pixel coordinates to celestial coordinates. The 
intermediate world coordinates of Paper I, Fig. 1 are here interpreted 
as projection plane coordinates, i.e. Cartesian coordinates in the plane 
of projection, and the multiple steps required to produce them have 
been condensed into one. This paper is concerned in particular with 
the steps enclosed in the dotted box. 



The original FITS paper by Wells et al. (1981) intro- 
duced the CRPIX ja l keyword to define the pixel coordinates 
( r i, f~2, r-}, ■ ■ .) of a coordinate reference point. Paper I retains 
this but replaces the coordinate rotation keywords CROTA/ with 
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1 The single-character alternate version code "a" on the various 
FITS keywords was introduced in Paper I. It has values blank and 
A through Z. 
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a linear transformation matrix. Thus, the transformation of 
pixel coordinates (p\,pi,p-i, ■ ■ ■) to intermediate world coor- 
dinates (x l ,x 2 , x 3 , . . .) becomes 

N 

Xi = v,^ »;,,(/); -r,). (1) 

7=1 

N 

7=1 

where N is the number of axes given by the NAX IS keyword. As 
suggested by the two forms of the equation, the scales s,-, and 
matrix elements m (i may be represented either separately or in 
combination. In the first form Sj is given by CDELT/a and m (J by 
PC i_ja; in the second, the product s,-m/y is given by CD i_ja. The 
two forms may not coexist in one coordinate representation. 

Equation (1) establishes that the reference point is the ori- 
gin of intermediate world coordinates. We require that the lin- 
ear transformation be constructed so that the plane of projec- 
tion is defined by two axes of the x, coordinate space. We will 
refer to intermediate world coordinates in this plane as pro- 
jection plane coordinates, (x,y), thus with reference point at 
(x,y) = (0, 0). Note that this does not necessarily correspond to 
any plane defined by the pj axes since the linear transformation 
matrix may introduce rotation and/or skew. 

Wells et al. (1981) established that all angles in FITS 
were to be measured in degrees and this has been entrenched 
by the AIPS convention and confirmed in the IAU-endorsed 
FITS standard (Hanisch et al., 2001). Paper I introduced the 
CUNIT/a keyword to define the units of CRVAL/a and CDELT/a. 
Accordingly, we require CUNIT/a = 'deg' for the celes- 
tial CRVAL/a and CDELT/a, whether given explicitly or not. 
Consequently, the (x, y) coordinates in the plane of projection 
are measured in degrees. For consistency, we use degree mea- 
sure for native and celestial spherical coordinates and for all 
other angular measures in this paper. 

For linear coordinate systems Wells et al. (1981) prescribed 
that world coordinates should be computed simply by adding 
the relative world coordinates, x,-, to the coordinate value at 
the reference point given by CRVAL ia. Paper I extends this by 
providing that particular values of CTYPE/a may be established 
by convention to denote non-linear transformations involving 
predefined functions of the x, parameterized by the CRVAL/a 
keyword values and possibly other parameters. 

In Sects. 2, 3 and 5 of this paper we will define the func- 
tions for the transformation from (x, y) coordinates in the plane 
of projection to celestial spherical coordinates for all spherical 
map projections likely to be of use in astronomy. 

The FITS header keywords discussed within the main body 
of this paper apply to the primary image header and image ex- 
tension headers. Image fragments within binary tables exten- 
sions defined by Cotton et al. (1995) have additional nomencla- 
ture requirements, a solution for which was proposed in Paper I. 
Coordinate descriptions may also be associated with the ran- 
dom groups data format defined by Greisen & Harten (1981). 
These issues will be expanded upon in Sect. 4. 

Section 6 considers the translation of older AIPS conven- 
tion FITS headers to the new system and provisions that may 



be made to support older FITS-reading programs. Section 7 
discusses the concepts presented here, including worked exam- 
ples of header interpretation and construction. 

2. Basic concepts 

2.1. Spherical projection 

As indicated in Fig. 1, the first step in transforming (x,y) co- 
ordinates in the plane of projection to celestial coordinates is 
to convert them to native longitude and latitude, (<p, 9). The 
equations for the transformation depend on the particular pro- 
jection and this will be specified via the CTYPE/a keyword. 
Paper I defined "4-3" form for such purposes; the rightmost 
three-characters are used as an algorithm code that in this pa- 
per will specify the projection. For example, the stereographic 
projection will be coded as 'STG'. Some projections require 
additional parameters that will be specified by the FITS key- 
words PV/_ma for m = 0, 1, 2, . . ., also introduced in Paper I. 
These parameters may be associated with the longitude and/or 
latitude coordinate as specified for each projection. However, 
definition of the three-letter codes for the projections and the 
equations, their inverses and the parameters which define them, 
form a large part of this work and will be discussed in Sect. 5. 
The leftmost four characters of CTYPE/a are used to identify 
the celestial coordinate system and will be discussed in Sect. 3. 

2.2. Reference point of the projection 

The last step in the chain of transformations shown in Fig. 1 is 
the spherical rotation from native coordinates, (<f>, 9), to celes- 
tial 2 coordinates (a, S). Since a spherical rotation is completely 
specified by three Euler angles it remains only to define them. 

In principle, specifying the celestial coordinates of any par- 
ticular native coordinate pair provides two of the Euler an- 
gles (either directly or indirectly). In the AIPS convention, the 
CRVAL ia keyword values for the celestial axes 3 specify the ce- 
lestial coordinates of the reference point and this in turn is as- 
sociated with a particular point on the projection. For zenithal 
projections that point is the sphere's point of tangency to the 
plane of projection and this is the pole of the native coordinate 
system. Thus the AIPS convention links a celestial coordinate 
pair to a native coordinate pair via the reference point. Note that 
this association via the reference point is purely conventional; 
it has benefits which are discussed in Sect. 5 but in principle 
any other point could have been chosen. 

Section 5 presents the projection equations for the transfor- 
mation of (x,y) to ((f), 9). The native coordinates of the refer- 
ence point would therefore be those obtained for (x, y) = (0, 0). 
However, it may happen that this point lies outside the bound- 
ary of the projection, for example as for the ZPN projection of 
Sect. 5.1.7. Therefore, while this work follows the AIPS ap- 
proach it must of necessity generalize it. 

2 Usage here of the conventional symbols for right ascension and 
declination for celestial coordinates is meant only as a mnemonic. It 
does not preclude other celestial systems. 

3 We will refer to these simply as "the CRVAL ia", and likewise for 
the other keyword values. 
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Accordingly we specify only that a fiducial celestial coor- 
dinate pair (ao,6o) given by the CRVAL/a will be associated 
with a fiducial native coordinate pair (fa, 9q) defined explicitly 
for each projection. For example, zenithal projections all have 
(0 O , O ) = (0, 90°), while cylindricals have (0 , 0o) = (0, 0). The 
AIPS convention has been honored here as far as practicable by 
constructing the projection equations so that (0o» #o) transforms 
to the reference point, (x,y) = (0, 0). Thus, apart from the one 
exception noted, the fiducial celestial and native coordinates 
are the celestial and native coordinates of the reference point 
and we will not normally draw a distinction. 

It is important to understand why ((/>q, 8o) differs for dif- 
ferent projection types. There are two main reasons; the first 
makes it difficult for it to be the same, while the second makes 
it desirable that it differs. Of the former, some projections such 
as Mercator's, diverge at the native pole, therefore they cannot 
have the reference point there because that would imply infinite 
values for CRPIX ja. On the other hand, the gnomonic projec- 
tion diverges at the equator so it can't have the reference point 
there for the same reason. Possibly (<p , 8 ) chosen at some mid- 
latitude could satisfy all projections, but that leads us to the 
second reason. 

Different projection types are best suited to different pur- 
poses. For example, zenithal projections are best for mapping 
the region in the vicinity of a point, often a pole; cylindrical 
projections are appropriate for the neighborhood of a great cir- 
cle, usually an equator; and the conies are suitable for small 
circles such as parallels of latitude. Thus, it would be awkward 
if a cylindrical used to map, say, a few degrees on either side 
of the galactic plane, had its reference point, and thus CRPIX ja 
and CRVAL ia, at the native pole, way outside the map bound- 
ary. In formulating the projection equations themselves the na- 
tive coordinate system is chosen to simplify the geometry as 
much as possible. For the zenithals the natural formulation has 
(x, y) = (0, 0) at the native pole, whereas for the cylindricals the 
equations are simplest if (x, y) = (0, 0) at a point on the equator. 

As discussed above, a third Euler angle must be specified 
and this will be given by the native longitude of the celestial 
pole, <f> p , specified by the new FITS keyword 

LONPOLEa (floating-valued). 

The default value of LONPOLEa will be for S Q > 6 or 180° 
for So < do- This is the condition for the celestial latitude to 
increase in the same direction as the native latitude at the refer- 
ence point. Thus, for example, in zenithal projections the de- 
fault is always 180° (unless So = 90°) since 8 Q = 90°. In 
cylindrical projections, where 6o = 0, the default value for 
LONPOLEa is for 6 Q > 0, but it is 180° for S < 0. 

2.3. Spherical coordinate rotation 

Since (<po, 9q) differs for different projections it is apparent that 
the relationship between (ao, So) and the required Euler angles 
also differs. 

For zenithal projections, (<po, 6o) = (0, 90°) so the CRVAL/a 
specify the celestial coordinates of the native pole, i.e. 
(ao,6o) = (a p ,S p ). There is a simple relationship between 



the Euler angles for consecutive rotations about the Z-, X- 
, and Z-axes and a p , S ? and (p v ; the ZXZ Euler angles are 
(a p + 90°, 90° -S p , (p p -90°). Given this close correspondence it 
is convenient to write the Euler angle transformation formulae 
directly in terms of a p , 6 P and (p p : 

a = a p + arg (sin cos £ p - cos sin <5 p cos(0 - tf> p ), 

- cos 8 sin(0 - <f> p )) , (2) 
S = sin _1 (sin0sin5 p + cos0cos5 p cos(0 - p )) . 

where arg () is an inverse tangent function that returns the cor- 
rect quadrant, i.e. if (x,y) = (rcos/3, r sin/3) with r > then 
arg(x,y) = /3. Note that, if 6 P = 90°, Eqs. (2) become 



a = ff p + 0-0p - 180°, 
6 = 8, 



(3) 



which may be used to define a simple change in the origin of 
longitude. Likewise for 6 P = -90° 



a = a v - <f> + <p p , 

5 = -e. 



(4) 



The inverse equations are 

(f> = (p p + arg (sin £ cos (5p - cos 5 sin £ p cos(a - a p ), 

- cos 6 sin(a - a p )) , (5) 
6 = sin _1 (sin5sin5 p + cos 5 cos 6 P cos(a - a p )) , 

Useful relations derived from Eqs. (2) and (5) are 

cos 6 cos(a - a p ) = sin 6 cos 6 P - cos 6 sin 6 P cos(0 - <p p ) , 

cos 6 sin(a - a p ) = - cos 6 sin(</> - (f> p ) , (6) 

cos 9cos((/> - <p p ) = sin 6 cos 6 P - cos 6 sin 6 P cos(a - a p ) , 
cos #sin(0 - <f> p ) = - cos 6 sin(a - a p ) . (7) 

A matrix method of handling the spherical coordinate rotation 
is described in Appendix A. 



2.4. Non-polar {(f>o,8o) 

Projections such as the cylindricals and conies for which 
((f>o, 0q) + (0, 90°) are handled by providing formulae to com- 
pute (or p , 5 p ) from (a , do) whence the above equations may be 
used. 

Given that (ao, do) are the celestial coordinates of the point 
with native coordinates (0o,#oX Eqs. (6) and (7) may be in- 
verted to obtain 

d p = arg (cos 8q cos(0 p - <po) , sin 8q) ± 



cos 



sin So 



■yjl - cos 2 8q sin 2 (0 p 



0o) 



sin(a'o - ap) = sin(0 p - <po) cos#o/ cos<5o , 
sin 8q - sin S p sin So 

cos(a -ap) = 



(8) 

(9) 
(10) 



cos Sp cos So 

whence Eqs. (2) may be used to determine the celestial coor- 
dinates. Note that Eq. (8) contains an ambiguity in the sign of 
the inverse cosine and that all three indicate that some com- 
binations of (po, 8o, ao, So, and <p v are not allowed. For these 
projections, we must therefore adopt additional conventions: 
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1. Equations (9) and (10) indicate that a p is undefined when 
do = ±90°. This simply represents the longitude singularity 
at the pole and forces us to define a p = ao in this case. 

2. If 6 p = ±90° then the longitude at the native pole is a p = 
a a + <A P ~ <Ao ~ 180° for 5 p = 90° and a p = a - p + <f> for 
6 P = -90°. 

3. Some combinations of <f>o, 6o, So, and p produce an in- 
valid argument for the cos _1 () in Eq. (8). This is indicative 
of an inconsistency for which there is no solution for S p . 
Otherwise Eq. (8) produces two solutions for S p . Valid so- 
lutions are ones that lie in the range -90° to +90°, and it is 
possible in some cases that neither solution is valid. 
Note, however, that if <f>o = 0, as is usual, then when 
LONPOLEa (= p ) takes its default value of or 180° (de- 
pending on 8q) then any combination of do and 9q produces 
a valid argument to the cos~'() in Eq. (8), and at least one 
of the solutions is valid. 

4. Where Eq. (8) has two valid solutions the one closest to the 
value of the new FITS keyword 

LATPOLEa (floating- valued). 

is chosen. It is acceptable to set LATPOLEa to a number 
greater than +90° to choose the northerly solution (the de- 
fault if LATPOLEa is omitted), or a number less than -90° 
to select the southern solution. 

5. Equation (8) often only has one valid solution (because 
the other lies outside the range -90° to +90°). In this case 
LATPOLEa is ignored. 

6. For the special case where 8q = 0, So = 0, and p - 0o = 
±90° then 5 P is not determined and LATPOLEa specifies it 
completely. LATPOLEa has no default value in this case. 

These rules governing the application of Eqs. (8), (9) and (10) 
are certainly the most complex of this formalism. FITS writers 
are well advised to check the values of <po, 6q, ao, So, and (p p 
against them to ensure their validity. 

2.5. User-specified (cf>o, %) 

In Sect. 2.2 we formally decoupled (ao, <5o) from the reference 
point and associated it with (<po, 6q). One implication of this is 
that it should be possible to allow (0o, # ) to be user-specifiable. 
This may be useful in some circumstances, mainly to allow 
CRVAL ia to match a point of interest rather than some prede- 
fined point which may lie well outside the image and be of no 
particular interest. We therefore reserve keywords PVLla and 
PV/_2a attached to longitude coordinate i to specify <f>o and 9q 
respectively. 

By itself, this prescription discards the AIPS convention 
and lacks utility because it breaks the connection between 
CRVAL ia and any point whose pixel coordinates are given in the 
FITS header. New keywords could be invented to define these 
pixel coordinates but this would introduce additional complex- 
ity and still not satisfy the AIPS convention. The solution 
adopted here is to provide an option to force (x, y) = (0, 0) at 
(<f>o, #o) by introducing an implied offset (xo,yo) which is com- 
puted for (0o, #o) from the relevant projection equations given 
in Sect. 5. This is to be applied to the (x,y) coordinates when 



converting to or from pixel coordinates. The operation is con- 
trolled by the value of PVL8a attached to longitude coordinate 
i; the offset is to be applied only when this differs from its de- 
fault value of zero. 

This construct should be considered advanced usage, of 
which Figs. 33 and 34 provide an example. Normally we ex- 
pect that PVLla and PV;'_2a will either be omitted or set to the 
projection-specific defaults given in Sect. 5. 

2.6. Encapsulation 

So that all required transformation parameters can be contained 
completely within the recognized world coordinate system 
(WCS) header cards, the values of LONPOLEa and LATPOLEa 
may be recorded as PVL3a and PVL4a, attached to longitude 
coordinate i, and these take precedence where a conflict arises. 

We recommend that FITS writers include the PVLla, 
PVL2a, PVL3a, and PVL4a cards in the header, even if only 
to denote the correct use of the default values. 

Note carefully that these are associated with the longitude 
coordinate, whereas the projection parameters defined later are 
all associated with the latitude coordinate. 

2.7. Change of coordinate system 

A change of coordinate system may be effected in a straightfor- 
ward way if the transformation from the original system, (a, S), 
to the new system, (a', S'), and its inverse are known. The new 
coordinates of the (0o> #o)> namely (a' , 6' ), are obtained simply 
by transforming (ao, So). To obtain <f>' first transform the coor- 
dinates of the pole of the new system to the original celestial 
system and then transform the result to native coordinates via 
Eq. (5) to obtain (<p' p , ff p ). As a check, compute 6' p via Eq. (8) 
and verify that 8' p = 5' p . 

2.8. Comparison with linear coordinate systems 

It must be stressed that the coordinate transformation described 
here differs from the linear transformation defined by Wells 
et al. (1981) even for some simple projections where at first 
glance they may appear to be the same. Consider the plate 
carree projection defined in Sect. 5.2.3 with = x, 9 = y and 
illustrated in Fig. 18. Figure 1 shows that while the transfor- 
mation from (x, y) to (0, 6) may be linear (in fact identical), 
there still remains the non-linear transformation from (0, 9) to 
(a, 6). Hence the linear coordinate description defined by the 
unqualified CTYPE/a pair of 'RA', 'DEC which uses the Wells 

et al. prescription will generally differ from that of 'RA CAR' , 

'DEC CAR' with the same CRVAL ia, etc. If LONPOLEa as- 
sumes its default value then they will agree to first order at 
points near the reference point but gradually diverge at points 
away from it. 

Figure 2 illustrates this point for a plate carree projection 
with reference coordinates of 8 hr right ascension and +60° dec- 
lination and with <p p = 0. It is evident that since the plate carree 
has (0o, do) = (0, 0), a non-oblique graticule may only be ob- 
tained by setting do = with <p p = 0. It should also be noted 
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Fig. 2. A linear equatorial coordinate system (top) defined via the 
methods of Wells et al. (1981), and the corresponding oblique system 
constructed using the methods of this paper. The reference coordinate 
(ao,do) for each is at right ascension 8 hr , declination +60° (marked). 
The two sets of FITS header cards differ only in their CTYPE ia key- 
word values. The non-oblique graticule could be obtained in the cur- 
rent system by setting (a , 6 Q , p ) = (120°, 0, 0). 

that where a larger map is to be composed of tiled submaps the 
coordinate description of a submap should only differ in the 
value of its reference pixel coordinate. 

3. Celestial coordinate systems 

As mentioned in Sect. 2.1, Paper I defined "4-3" form for the 
CTYPE ia keyword value; i.e., the first four characters specify 
the coordinate type, the fifth character is a and the remain- 
ing three characters specify an algorithm code for computing 
the world coordinate value. Thus while the right half specifies 
the transformation to be applied in computing the spherical co- 
ordinate pair, (a, 6), the left half simply identifies the celes- 
tial system with which (a, 6) are to be associated. In this sense 
CTYPE ia contains an active part which drives the transforma- 
tion and a passive part which describes the results. 

Consistent with past practice, an equatorial coordinate pair 
is denoted by 'RA--' and 'DEC-', and other celestial systems 
are of the form 'xLON' and 'xLAT' for longitude and latitude 
pairs, where x = G for galactic 4 , E for ecliptic, H for helioe- 
cliptic 5 , and S for supergalactic coordinates. Since represen- 

4 "New" galactic coordinates are assumed here, Blaauw et al. 
(1960). Users of the older system or future systems should adopt a 
different value of x and document its meaning. 

5 Ecliptic and helioecliptic systems each have their equator on the 
ecliptic. However, the reference point for ecliptic longitude is the ver- 
nal equinox while that for helioecliptic longitude is the sun vector. 



tation of planetary, lunar and solar coordinate systems could 
exceed the 26 possibilities afforded by the single character x, 
we also here allow 'yzLN' and 'yzLT pairs. Additional values 
of x and yz will undoubtedly be defined. A basic requirement, 
however, is that the coordinate system be right-handed and that 
the pole be at latitude +90°. A coordinate pair such as azimuth 
and zenith distance would have to be represented as a negative 
azimuth and elevation with implied conversion. In some situa- 
tions negation of the longitude coordinate may be implemented 
via a sign reversal of the appropriate CDELT/a. It remains the 
responsibility of the authors of new coordinate system types to 
define them properly and to gain general recognition for them 
from the FITS community. However, FITS interpreters should 
be able to associate coordinate pairs even if the particular coor- 
dinate system is not recognized. 

Paper I clarified that, while the subscript of CRPIX ja and 
the j subscript of PC ija and CD i_ja refer to pixel coordinate el- 
ements, the i subscript of PCLy'a and CDLja and the subscripts 
of CDELT/a, CTYPE ia and CRVAL/a refer to world coordinate 
elements. However we now have three different sets of world 
coordinates, (x,y), {<p,ff), and (a, 6). This leads us to associate 
x, 0, and a on the one hand, and y, 9, and 6 on the other. This 
simply means, for example, that if CTYPE3 = 'GL0N-AIT', 
then the third element of the intermediate world coordinate cal- 
culated via Eq. (1) corresponds to what we have been calling 
the x-coordinate in the plane of projection, the association be- 
ing between a and x. In this way pairs of CTYPE ia with com- 
plementary left halves and matching right halves define which 
elements of the intermediate world coordinate vector form the 
plane of projection. 

3.1. Equatorial and ecliptic coordinates 

Several systems of equatorial coordinates (right ascension and 
declination) are in common use. Apart from the International 
Celestial Reference System (ICRS, IAU 1998), the axes of 
which are by definition fixed with respect to the celestial 
sphere, each system is parameterized by time. In particular, 
mean equatorial coordinates (Hohenkerk et al., 1992) are de- 
fined in terms of the epoch (i.e. instant of time) of the mean 
equator and equinox (i.e. pole and origin of right ascension). 
The same applies for ecliptic coordinate systems. Several addi- 
tional keywords are therefore required to specify these systems 
fully. We introduce the new keyword 

RADESYSa (character- valued) 

to specify the particular system. Recognized values are given in 
Table 1. Apart from FK4-N0-E these keywords are applicable 
to ecliptic as well as equatorial coordinates. 

Wells et al. (1981) introduced the keyword EPOCH to mean 
the epoch of the mean equator and equinox. However we here 
replace it with 

EQUINOX a (floating-valued), 
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Table 1. Allowed values of RADESYSa. 



ivfiUHij I ju 




'ICRS' 


International Celestial Reference System 


'FK5' 


mean place, 




new (IAU 1984) system 


'FK4' 


mean place, 




old (Bessell-Newcomb) system 


'FK4-N0-E' 


mean place, 




old system but without e-terms 


'GAPPT' 


geocentric apparent place, 




IAU 1984 system 



since the word "epoch" is often used in astrometry to refer to 
the time of observation. The new keyword 6 takes preference 
over EPOCH if both are given. Note that EQUINOXa applies to 
ecliptic as well as to equatorial coordinates. 

For RADESYSa values of FK4 and FK4-N0-E, any stated 
equinox is Besselian and, if neither EQUINOXa nor EPOCH are 
given, a default of 1950.0 is to be taken. For FK5, any stated 
equinox is Julian and, if neither keyword is given, it defaults to 
2000.0. 

If the EQUINOXa keyword is given it should always be ac- 
companied by RADESYSa. However, if it should happen to ap- 
pear by itself then RADESYSa defaults to FK4 if EQUINOXa 
< 1984.0, or to FK5 if EQUINOXa > 1984.0. Note that these 
defaults, while probably true of older files using the EPOCH key- 
word, are not required of them. 
RADESYSa 

defaults to ICRS if both it and EQUINOXa are absent. 

Geocentric apparent equatorial and ecliptic coordinates 
(GAPPT) require the epoch of the equator and equinox of date. 
This will be taken as the time of observation rather than 
EQUINOXa. The time of observation may also be required for 
other astrometric purposes in addition to the usual astrophys- 
ical uses, for example, to specify when the mean place was 
correct in accounting for proper motion, including "fictitious" 
proper motions in the conversion between the FK4 and FK5 sys- 
tems. The old DATE-OBS keyword may be used for this pur- 
pose. However, to provide a more convenient specification we 
here introduce the new keyword 

MJD-OBS (floating-valued), 

that provides DATE-OBS as a Modified Julian Date (JD - 
2400000.5) but is identical to it in all other respects. MJD-OBS 
does not have a version code since there can only be one time 
of observation. Following the year-2000 conventions for DATE 
keywords (Bunclark & Rots 1996), this time refers by default 
to the start of the observation unless another interpretation is 
clearly explained in the comment field. In the present case the 
distinction is not important. We leave it to future agreements to 
clarify systems of time measurement and other matters related 
to time. 

6 EQUINOX has since been adopted by Hanisch et al. (2001) which 
deprecates EPOCH. 



Table 2. Alternate coordinate keywords; the data type of the alternate 
keyword matches that of the primary keyword 



Keyword 


Primary 


BINTABLE 


Pixel 


Description 


Array 


Array 


List 


Coord Rotation 


LONPOLEa 


LONPna 


LONPna 


Coord Rotation 


LATPOLEa 


LATP na 


LATP na 


Coord Epoch 


EQUINOXa 


EQUIna 


EQUIna 


Date of Obs 


MJD-OBS 


MJDOBrc 


MJDOBrc 


Reference Frame 


RADESYSa 


RADEna 


RADEna 



The combination of CTYPE ia, RADESYSa, and EQUINOXa 
define the coordinate system of the CRVAL/a and of the ce- 
lestial coordinates that result from the sequence of transfor- 
mations summarized by Fig. 1. However, FK4 coordinates are 
not strictly spherical since they include a contribution from 
the elliptic terms of aberration, the so-called "e-terms" which 
amount to < 343 milliarcsec. Strictly speaking, therefore, a 
map obtained from, say, a radio synthesis telescope, should be 
regarded as FK4-N0-E unless it has been appropriately resam- 
pled or a distortion correction provided (Paper IV). In common 
usage, however, CRVAL ia for such maps is usually given in FK4 
coordinates. In doing so, the e-terms are effectively corrected 
to first order only. 

4. Alternate FITS image representations 

4.1. Random groups visibility data 

The random-groups extension to FITS (Greisen & Harten 
1981) has been used to transmit interferometer fringe visibility 
sample data. It has been customary among users of this format 
to convey the suggested projection type as the last four char- 
acters of the random parameter types (PTYPEn) of the Fourier 
plane coordinates (u, v, w). Any rotation of these coordinates 
was carried, however, by CR0TA; associated with the one -pixel 
array axis used to convey the field declination. Within the new 
convention users of random groups for visibility data should 
be prepared to read, carry forward, and use as needed, the 
new keywords PCLja, CDLja, LONPOLEa, LATPOLEa, PVLja, 
MJD-OBS, EQUINOXa, and RADESYSa. In particular, any rota- 
tion of the (u, v) coordinate system should be represented via 
PC i_ja or CD Lja (but not both) attached to the two degenerate 
axes used to convey the celestial coordinates of the field center. 

4.2. Pixel list and image array column elements 

Paper I defined the coordinate keywords used to describe a 
multi-dimensional image array in a single element of a FITS 
binary table and a tabulated list of pixel coordinates in a FITS 
ASCII or binary table. In this section we extend this to the key- 
words specific to celestial coordinates. 

4.2.1. Keyword naming convention 

Table 2 lists the corresponding set of coordinate system key- 
words for use with each type of FITS image representation. 
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The data type of the alternate keyword matches that of the 
primary keyword and the allowed values are the same. The fol- 
lowing notes apply to the naming conventions used in Table 2: 

- a is a 1 -character coordinate version code and may be blank 
(for the primary version of the coordinate description) or 
any single uppercase character from A through Z. 

- n is an integer table column number without any leading 
zeros (1 -999). 

When using the BINTABLE image array format (Cotton et 
al., 1995), if the table only contains a single image column or 
if there are multiple image columns but they all have the same 
value for any of the keywords in Table 2 then the simpler form 
of the keyword name, as used for primary arrays, may be used. 
For example, if all the images in the table have the same epoch 
then one may use a single EQUINOX a keyword rather than mul- 
tiple EQUI na keywords. The other keywords, however, must al- 
ways be specified using the more complex keyword name with 
the column number suffix and the axis number prefix. 

In principle, more than one pixel list image can be stored 
in a single FITS table by defining more than one pair of p\ 
and p2 pixel coordinate columns. Under the convention defined 
here, however, all the images must share the same values for the 
keywords in Table 2. 

Example binary table and pixel list headers are given in 
Sect. 7.3.2 

5. Spherical map projections 

In this section we present the transformation equations for all 
spherical map projections likely to be of use in astronomy. 
Many of these such as the gnomonic, orthographic, zenithal 
equidistant, Sanson-Flamsteed, Hammer-Aitoff and COBE 
quadrilateralized spherical cube are in common use. Others 
with special properties such as the stereographic, Mercator, and 
the various equal area projections could not be excluded. A se- 
lection of the conic and polyconic projections, much favored by 
cartographers for their minimal distortion properties, has also 
been included. However, we have omitted numerous other pro- 
jections which we considered of mathematical interest only. 
Evenden (1991) presents maps of the Earth for 73 different 
projections, although without mathematical definition, includ- 
ing most of those described here. These are particularly useful 
in judging the distortion introduced by the various projections. 
Snyder (1993) provides fascinating background material on the 
subject; historical footnotes in this paper, mainly highlighting 
astronomical connections, are generally taken from this source. 
It should be evident from the wide variety of projections de- 
scribed here that new projections could readily be accommo- 
dated, the main difficulty being in obtaining general recogni- 
tion for them from the FITS community. 

Cartographers have often given different names to spe- 
cial cases of a class of projection. This applies particularly to 
oblique projections which, as we have seen in Sect. 2, the cur- 
rent formalism handles in a general way. While we have tended 
to avoid such special cases, the gnomonic, stereographic, and 
orthographic projections, being specializations of the zenithal 



perspective projection, are included for conformance with the 
AIPS convention. It is also true that zenithal and cylindrical 
projections may be thought of as special cases of conic projec- 
tions (see Sect. 5.4). However, the limiting forms of the conic 
equations tend to become intractable and infinite-valued pro- 
jection parameters may be involved. Even when the conic equa- 
tions don't have singularities in these limits it is still likely to 
be less efficient to use them than the simpler special-case equa- 
tions. Moreover, we felt that it would be unwise to disguise 
the true nature of simple projections by implementing them as 
special cases of more general ones. In the same vein, the cylin- 
drical equal area projection, being a specialization of the cylin- 
drical perspective projection, stands on its own right, as does 
the Sanson-Flamsteed projection which is a limiting case of 
Bonne's projection. A list of aliases is provided in Appendix C, 
Table A. 1. 

The choice of a projection often depends on particular spe- 
cial properties that it may have. Certain equal area projections 
(also known as authalic, equiareal, equivalent, homalographic, 
homolo graphic, or homeotheric) have the property that equal 
areas on the sphere are projected as equal areas in the plane 
of projection. This is obviously a useful property when surface 
density must be preserved. Conformality is a property which 
applies to points in the plane of projection which are locally 
distortion-free. Practically speaking, this means that the pro- 
jected meridian and parallel through the point intersect at right 
angles and are equiscaled. A projection is said to be confor- 
mal or orthomorphic if it has this property at all points. Such 
a projection cannot be equal area. It must be stressed that con- 
formality is a local property, finite regions in conformal projec- 
tions may be very severely distorted. A number of projections 
have other special properties and these will be noted for each. 

Maps of the Earth are conventionally displayed with ter- 
restrial latitude increasing upwards and longitude to the right, 
i.e. north up and east to the right, as befits a sphere seen from 
the outside. On the other hand, since the celestial sphere is 
seen from the inside, north is conventionally up and east to the 
left. The AIPS convention arranged that celestial coordinates at 
points near the reference point should be calculable to first or- 
der via the original linear prescription of Wells et al. (1981), i.e. 
(a, 6) « (ao,<5o) + (x,y). Consequently, the CDELT/a keyword 
value associated with the right ascension was negative while 
that for the declination was positive. The handedness of the 
(x, y) coordinates as calculated by the AIPS convention equiv- 
alent of Eq. ( 1 ) is therefore opposite to that of the (p i , pi) pixel 
coordinates. 

In accordance with the image display convention of Paper I 
we think of the p\ -pixel coordinate increasing to the right with 
P2 increasing upwards, i.e. a right-handed system. This means 
that the (x,y) coordinates must be left-handed as shown in 
Fig. 3. Note, however, that the approximation (a, 6) ~ (ao, Sq) + 
(x,y) cannot hold unless 1) (<f>o,6o), and hence (ao,So), d° ac ~ 
tually map to the reference point (Sect. 2.2), 2) (f> p assumes its 
default value (Sect. 2.8), and 3) the projection is scaled true 
at the reference point (some are not as discussed in Sect. 7.2). 
Figure 3 also illustrates the orientation of the native coordinate 
system with respect to the (x, y) coordinate system for the two 
main cases. 
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Cartographers, for example Kellaway (1946), think of 
spherical projections as being a projection of the surface of a 
sphere onto a plane, this being the forward direction; the depro- 
jection from plane back to sphere is thus the inverse or reverse 
direction. However, this is at variance with common usage in 
FITS where the transformation from pixel coordinates to world 
coordinates is considered the forward direction. We take the 
cartographic view in this section as being the natural one and 
trust that any potential ambiguity may readily be resolved by 
context. 

The requirement stated in Sect. 1 that (x, y) coordinates in 
the plane of projection be measured in "degrees" begs clari- 
fication. Spherical projections are usually defined mathemati- 
cally in terms of a scale factor, ro, known as the "radius of the 
generating sphere". However, in this work ro is set explicitly 
to 180°/7r in order to maintain backwards compatibility with 
the AIPS convention. This effectively sets the circumference of 
the generating sphere to 360° so that arc length is measured 
naturally in degrees (rather than radians as for a unit sphere). 
However, this true angular measure on the generating sphere 
becomes distorted when the sphere is projected onto the plane 
of projection. So while the "degree" units of ro are notionally 
carried over by conventional dimensional analysis to the (x,y) 
they no longer represent a true angle except near the reference 
point (for most projections). 

In addition to the (x, y) coordinates, the native spherical co- 
ordinates, ((f), ff), celestial coordinates, (a, 6), and all other an- 
gles in this paper are measured in degrees. In the equations 
given below, the arguments to all trigonometric functions are 
in degrees and all inverse trigonometric functions return their 
result in degrees. Whenever a conversion between radians and 
degrees is required it is shown explicitly. All of the graticules 
presented in this section have been drawn to the same scale 
in (x, y) coordinates in order to represent accurately the exag- 
geration and foreshortening found in some projections. It will 
also be apparent that since FITS image planes are rectangular 
and the boundaries of many projections are curved, there may 
sometimes be cases when the FITS image must contain pixels 



that are outside the boundary of the projection. These pixels 
should be blanked correctly and geometry routines should re- 
turn a sensible error code to indicate that their celestial coordi- 
nates are undefined. 

5.1. Zenithal (azimuthal) projections 

Zenithal or azimuthal projections all map the sphere directly 
onto a plane. The native coordinate system is chosen to have 
the polar axis orthogonal to the plane of projection at the refer- 
ence point as shown on the left side of Fig. 3. Meridians of na- 
tive longitude are projected as uniformly spaced rays emanat- 
ing from the reference point and the parallels of native latitude 
are mapped as concentric circles centered on the same point. 
Since all zenithal projections are constructed with the pole of 
the native coordinate system at the reference point we set 

(00, ^zenithal = (0,90°). (11) 

Zenithal projections are completely specified by defining 
the radius as a function of native latitude, Rg. Rectangular 
Cartesian coordinates, (x,y), in the plane of projection as de- 
fined by Eq. (1), are then given by 

x = Rgsin(p, (12) 
y = -Rgcoscp, (13) 

which may be inverted as 

<p = arg(-y,x), (14) 

R e = ^+y 2 . (15) 

5.1 .1 . AZP: zenithal perspective 

Zenithal (azimuthal) perspective projections are generated 
from a point and carried through the sphere to the plane of 
projection as illustrated in Fig. 4. We need consider only the 
case where the plane of projection is tangent to the sphere at 
its pole; the projection is simply rescaled if the plane intersects 
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some other parallel of native latitude. If the source of the pro- 
jection is at a distance /j. spherical radii from the center of the 
sphere with /j. increasing in the direction away from the plane 
of projection, then it is straightforward to show that 



Re 



180° (yu+ l)cos0 
n [i + sin 9 



(16) 



with /i + -1 being the only restriction. When Rg is given 
Eq. (16) has two solutions for 9, one for each side of the sphere. 
The following form of the inverse equation always gives the 
planeward solution for any 



9 = arg(p, 1) - sin 



p/j. 



where 



Rn 



P = 



180° fi+ 1 



(17) 



(18) 



For |yu| + 1 the sphere is divided by a native parallel at latitude 
9 X into two unequal segments that are projected in superposi- 
tion: 



9, 



sm\-l/fx) ...\fx\> 1 
sin _1 (-/0 ...\fx\<l 



(19) 



For \fx\ > 1, the projection is bounded and both segments are 
projected in the correct orientation, whereas for < 1 the 
projection is unbounded and the anti-planewards segment is in- 
verted. 

A near-sided perspective projection may be obtained with 
fi < -1. This correctly represents the image of a sphere, such 
as a planet, when viewed from a distance times the planetary 
radius. The coordinates of the reference point may be expressed 
in planetary longitude and latitude, (A.,/3). Also, the signs of the 
relevant CDELT/a may be chosen so that longitude increases as 
appropriate for a sphere seen from the outside rather than from 
within. 

It is particularly with regard to planetary mapping that we 
must generalize AZP to the case where the plane of projection is 
tilted with respect to the axis of the generating sphere, as shown 
on the left side of Fig. 5. It can be shown (Sect. 7.4.1) that this 
geometry is appropriate for spacecraft imaging with non-zero 



look-angle, y, the angle between the camera's optical axis and 
the line to the center of the planet. 

Such slant zenithal perspective projections are not radially 
symmetric and their projection equations must be expressed di- 
rectly in terms of x and y: 



x — Rsincfi, 

y = -R secy cos cp , 

where 

180° 0u + l)cos6» 



R = 



(20) 
(21) 

— (22) 
n (p, + sin 9) + cos 9 cos cp tan y 

is a function of <f> as well as 9. Equations (20) and (21) reduce 
to the non-slant case for y — 0. The inverse equations are 

cp = arg(-ycosy, x) , (23) 

iff - 0) 

iA + w + 180 ' 

where 



(24) 



CO 



arg(p, 1), 
( 

sin -1 



2 + l 

R 



P = 



R 



if-(^ + l)+ysiny 



jx 2 + y 2 cos 2 y. 



(25) 
(26) 

(27) 
(28) 



For \fi\ < 1 only one of the solutions for 9 will be valid, i.e. lie in 
the range [-90°, 90°] after normalization. Otherwise there will 
be two valid solutions; the one closest to 90° should be chosen. 

With y + the projection is not scaled true at the reference 
point. In fact the x scale is correct but the y scale is magnified 
by sec y, thus stretching parallels of latitude near the pole into 
ellipses (see Fig. 6). This also shows the native meridians pro- 
jected as rays emanating from the pole. For constant 9, each 
parallel of native latitude defines a cone with apex at the point 
of projection. This cone intersects the tilted plane of projec- 
tion in a conic section. Equations (20) and (21) reduce to the 
parametric equations of an ellipse, parabola, or hyperbola; the 
quantity 

C = (p + sin 9) 2 - tan 2 y cos 2 9 (29) 
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Fig. 5. Alternate geometries of slant zenithal perspective projections with ji = 2 and y = 30°: (left) tilted plane of projection (AZP), (right) 
displaced point of projection P (SZP). Grey lines in each diagram indicate the other point of view. They differ geometrically only by a scale 
factor, the effect of translating the plane of projection. Each projection has its native pole at the reference point, r and r', but these are 
geometrically different points. Thus the native latitudes 6 and ff of the geometrically equivalent points s and s' differ. Consequently the two sets 
of projection equations have a different form. This diagram is at 3:2 scale compared to the graticules of Figs. 6 and 7. 



determines which: 
C > ...ellipse, 

C = ...parabola, (30) 
C < ...hyperbola. 

If |/i cos y\ < 1 then the open conic sections are possible and 
C = when 



±y - sin '(/icosy) . 



(31) 



For C > the eccentricity of the ellipse is a function of 9, as is 
the offset of its center in y. 

Definition of the perimeter of the projection is more com- 
plicated for the slant projection than the orthogonal case. As 
before, for \fi\ > 1 the sphere is divided into two unequal seg- 
ments that are projected in superposition. The boundary be- 
tween these two segments is what would be seen as the limb of 
the planet in spacecraft photography. It corresponds to native 
latitude 



Ax 



sin _1 (- 



(32) 



which is projected as an ellipse, parabola, or hyperbola for 
I// cos y\ greater than, equal to, or less than 1 respectively. 

In general, for \ft cos y\ > 1, the projection is bounded, oth- 
erwise it is unbounded. However, the latitude of divergence is 
now a function of <p: 



where 



if/ - h) 

tf/ + oj+ 180° 



L) = 



■tan ^tanycos^), 



sin 



(33) 

(34) 
(35) 



, Vl + tan 2 y cos 2 ip 

Zero, one, or both of the values of given by Eq. (33) may be 
valid, i.e. lie in the range [-90°, 90°] after normalization. 

The FITS keywords PVLla and PV/_2a, attached to latitude 
coordinate i, will be used to specify, respectively, [i in spherical 
radii and y in degrees, both with default value 0. 



5.1 .2. SZP: slant zenithal perspective 

While the generalization of the AZP projection to tilted planes 
of projection is useful for certain applications it does have a 
number of drawbacks, in particular, unequal scaling at the ref- 
erence point. 

Figure 5 shows that moving the point of projection, P, off 
the axis of the generating sphere is equivalent, to within a scale 
factor, to tilting the plane of projection. However this approach 
has the advantage that the plane of projection remains tangent 
to the sphere. Thus the projection is conformal at the native 
pole as can be seen by the circle around the native pole in Fig. 7. 
It is also quite straightforward to formulate the projection equa- 
tions with P offset in x as well as y. 

It is interesting to note that this slant zenithal perspec- 
tive (SZP) projection also handles the case that corresponds to 
y = 90° in AZP. AZP fails in this extreme since P falls in the 
plane of projection - effectively a scale factor of zero is ap- 
plied to AZP over the corresponding SZP case. One of the more 
important aspects of SZP is the application of its limiting case 
with ji — co in aperture synthesis radio astronomy as discussed 
in Sect. 5. 1 .5. One minor disadvantage is that the native merid- 
ians are projected as curved conic sections rather than straight 
lines. 

If the Cartesian coordinates of P measured in units of the 
spherical radius are (x p ,y p ,z p ), then 

1 80° z P cos 6 sin <p - x p ( 1 - sin 6) 



n z p - (1 - sin#) 

1 80° z P cos cos <f> + y p ( 1 - sin ff) 
n z p - ( 1 - sin 0) 

To invert these equations, compute 6 first via 



6 



sin 



where 



-b± VP" 



X' 2 + Y' 2 + 1 , 



(36) 



(37) 



(38) 



(39) 
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Fig. 6. Slant zenithal perspective (AZP) projection with li = 2 and 
7 = 30° which corresponds to the left-hand side of Fig. 5. The 
reference point of the corresponding SZP projection is marked at 
(4>,0) = (0,60°). 



081 








Fig. 7. Slant zenithal perspective (SZP) projection with fi = 2 and 
(4> C ,8 C ) = (180°, 60°), which corresponds to the right-hand side of 
Fig. 5. The reference point of the corresponding AZP projection is 
marked at (<p, 6) = (180°, 60°). 



b = X'(X -X') + Y'(Y -Y'), 
c = (X - X'f + (Y- Y'f - 1 , 



(X,Y) 



180' 



:(x,y), 



(X',Y') = (X-x p ,Y-y p )/ Zp . 



(40) 
(41) 
(42) 
(43) 



Choose 6 closer to 90° if Eq. (38) has two valid solutions; then 
(p is given by 



tp = arg (-(Y - y'(l - sin 6)), X-X'(l- sin 0)) . 



(44) 



The Cartesian coordinates of P are simply related to the pa- 
rameters used in AZP. If li is the distance of P from the center of 
the sphere O, and the line through P and O intersects the sphere 
at (<f> c , 6 C ) on the planewards side (point r in Fig. 5, right), then 
y = 90° - 6 C and 



x p = -ft cos # c sin <f> c , 
y ? = pL cos 6 C cos cb c , 
z p = li sin 6> c + 1 , 



(45) 
(46) 
(47) 



where li is positive if P lies on the opposite side of O from 
r and negative otherwise. For a non-degenerate projection we 
require z p + and this is the only constraint on the projection 
parameters. 

For > 1 the sphere is divided into two unequal segments 
that are projected in superposition. The limb is defined by com- 
puting the native latitude 9 X as a function of $ 



if/ — CO 

\p + oj+ 180 ' 
where 

iff = arg(p,cr), 

co = sin" 



-y/p 2 +cr 2 j 



(48) 

(49) 
(50) 



(p,cr) = (z p - 1, x p sin0 -y p cos<f>), 

= iji sin 6 C , -j-i cos 8 C cos(0 - cb c ) . 



(51) 
(52) 



Zero, one, or both of the values of 9 X given by Eq. (48) may 
be valid, i.e. lie in the range [-90°, 90°] after normalization. A 
second boundary constraint applies if 1 1 - z p | < 1 in which case 
the projection diverges at native latitude: 



Boo = sin '(1 -z p ). 



(53) 



The FITS keywords PVLla, PV/_2a, and PV/_3a, attached 
to latitude coordinate i, will be used to specify, respectively, li 
in spherical radii with default value 0, cp c in degrees with default 
value 0, and 9 C in degrees with default value 90°. 

5.1.3. TAN: gnomonic 

The zenithal perspective projection with li — 0, the gnomonic 
projection 7 , is widely used in optical astronomy and was given 
its own code within the AIPS convention, namely TAN 8 . For 
jj. = 0, Eq. (16) reduces to 



180° 
Re = cot 6 , 



with inverse 



6 = tan 



180' 



(54) 



(55) 



The gnomonic projection is illustrated in Fig. 8. Since the pro- 
jection is from the center of the sphere, all great circles are 



7 The gnomonic projection is the oldest known, dating to Thales of 
Miletus (ca. 624-547 B.C.). The stereographic and orthographic date to 
Hipparchus (ca. 190-after 126 B.C.). 

8 Referring to the dependence of R e on the angular separation be- 
tween the tangent point and field point, i.e. the native co-latitude. 
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Fig. 8. Gnomonic (TAN) projection; diverges at 6 = 0. 




■30 V> v 4/ ^ 

Fig. 9. Stereographic (STG) projection; diverges at 6 = -90°. 

projected as straight lines. Thus, the shortest distance between 
two points on the sphere is represented as a straight line inter- 
val, which, however, is not uniformly divided. The gnomonic 
projection diverges at 6 = 0, but one may use a gnomonic pro- 
jection onto the six faces of a cube to display the whole sky. 
See Sect. 5.6.1 for details. 



5.1.4. STG: stereographic 

The stereographic projection, the second important special case 
of the zenithal perspective projection defined by the AIPS con- 
vention, has ji = 1, for which Eq. (16) becomes 



Ra = 



180° 2cos0 
n 1 + sin 8 




Fig. 10. Slant orthographic (SIN) projection: (top) the orthographic 
projection, (£, if) = (0, 0), north and south sides begin to overlap at 6 = 
0; (bottom left) (0 C ,0 C ) = (225°, 60°), i.e. (£17) = (-1/ V6, 1/ V6); 
(bottom right) projection appropriate for an east-west array observing 
at So = 60°, (A, 0c) = (180°, 60°), (£ r,) = (0, 1/ V5). 



360° /90 o - 6 
tan 



with inverse 



= 90° -2 tan 



1 360°/ 



(57) 



The stereographic projection illustrated in Fig. 9 is the confor- 
mal (orthomorphic) zenithal projection 9 , everywhere satisfying 
the isoscaling requirement 

dRe = -nRg 
86 180°cos6>' 



(58) 



It also has the amazing property that it maps all circles on the 
sphere to circles in the plane of projection, although concentric 
circles on the sphere are not necessarily concentric in the plane 
of projection. This property made it the projection of choice for 
Arab astronomers in constructing astrolabes. 

5.1.5. SIN: slant orthographic 

The zenithal perspective projection with fj. = 00, the ortho- 
graphic projection, is illustrated in the upper portion of Fig. 10 
(at consistent scale). It represents the visual appearance of a 
sphere, e.g. a planet, when seen from a great distance. 

The orthographic projection is widely used in aperture syn- 
thesis radio astronomy and was given its own code within the 
AIPS convention, namely SIN 10 . Use of this projection code 
obviates the need to specify an infinite value as a parameter of 
AZP. In this case, Eq. (16) becomes 



R* = 



180° 



n 



cos 6 , 



(59) 



(56) 



9 First noted by astronomer Edmond Halley (1656-1742). 

10 Similar etymology to TAN. 
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Fig. 11. Zenithal equidistant (ARC) projection; no limits. 



Fig. 12. Zenithal polynomial projection (ZPN) with parameters, 0.050, 
0.975, -0.807, 0.337, -0.065, 0.010, 0.003, -0.001; limits depend upon 
the parameters. 



with inverse 

(l^ 9 ) 



6 = cos 



(60) 



In fact, use of the orthographic projection in radio inter- 
ferometry is an approximation, applicable only for small field 
sizes. However, an exact solution exists where the interferome- 
ter baselines are co-planar. It reduces to what Greisen (1983) 
called the NCP projection for the particular case of an East- 
West interferometer (Brouw 1971). The projection equations 
(derived in Appendix B) are 

180° r 

■ [cos#sin0 + £(1 - sin#)J , 



x = 



y = — 



n 
180° 



[cos#cos0 - rj(l - sin#)] 



(61) 
(62) 



These are the equations of the "slant orthographic" projection, 
equivalent to Eqs. (36) and (37) of the SZP projection in the 
limit fi — co, with 



n 



cot 6 C sin (f> c , 
• cot 6 C cos C . 



(63) 
(64) 



It can be shown that the slant orthographic projection is 
equivalent to an orthographic projection centered at (x,y) = 
rj) which has been stretched in the C direction by a fac- 
tor of cosec 9 C . The projection equations may be inverted using 
Eqs. (38) and (44) except that Eq. (43) is replaced with 

(X',r) = (£/?). (65) 

The outer boundary of the SIN projection is given by Eq. (48) 
in the limit fi = oo: 

8 X = - tan -1 (£ sin <p - rj cos <p) . (66) 



Two example graticules are illustrated in the lower portion 
of Fig. 10. We here extend the original SIN projection of the 
AIPS convention to encompass the slant orthographic projec- 
tion, with the dimensionless £ and rj given by keywords PVLla 
and PV;'_2a, respectively, attached to latitude coordinate i, both 
with default value 0. 



5.1.6. ARC: zenithal equidistant 

Some non-perspective zenithal projections are also of interest 
in astronomy. The zenithal equidistant projection first appeared 
in Greisen (1983) as ARC. It is widely used as the approximate 
projection of Schmidt telescopes. As illustrated in Fig. 11, the 
native meridians are uniformly divided to give equispaced par- 
allels. Thus 



R e = 90° -6, 



(67) 



which is trivially invertible. This projection was also known in 
antiquity. 

5.1.7. ZPN: zenithal polynomial 

The zenithal polynomial projection, ZPN, generalizes the ARC 
projection by adding polynomial terms up to a large degree in 
the zenith distance. We define it as 



m=0 



(68) 



Note the dimensionless units of P m imparted by 7r/180°. 

Since its inverse cannot be expressed analytically, ZPN 
should only be used when the geometry of the observations 
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Fig. 13. Zenithal equal area projection (ZEA); no limits. 



require it. In particular, it should never be used as an n th -degree 
expansion of one of the standard zenithal projections. 

If Pq is non-zero the native pole is mapped to an open circle 
centered on the reference point as illustrated in Fig. 12. In other 
words, (0 O , 6 ) = (0, 90°) is not at (x,y) = (0, 0), which in fact 
lies outside the boundary of the projection. However, we do 
not dismiss Po + as a possibility since it is not inconsistent 
with the formalism presented in Sect. 2.2 and could conceiv- 
ably be useful for images which do not contain the reference 
point. Needless to say, care should be exercised in constructing 
and interpreting such systems particularly in that (ao, So) (i.e. 
the CRVAL id) do not specify the celestial coordinates of the ref- 
erence point (the CRPIX ja). 

P m (dimensionless) is given by the keywords PV/JSa, 
PV/_la, . . ., PV/_29a, attached to latitude coordinate i, all of 
which have default values of zero. 



5.1.8. ZEA: zenithal equal-area 

Lambert's zenithal equal-area projection illustrated in Fig. 13 
is constructed by defining R g so that the area enclosed by the 
native parallel at latitude 8 in the plane of projection is equal to 
the area of the corresponding spherical cap. It may be generated 
using 



Re 



180° 

n 
360° 



V2(l -sin0) 



(69) 



sin 



90° -0 



with inverse 



= 90° -2 sin" 



1 360°/ ' 



(70) 




Fig. 14. Airy projection (AIR) with 6 h = 45°; diverges at 6 = -90°. 

5.1.9. AIR: Airy projection 

The Airy projection 1 1 minimizes the error for the region within 
latitude 8b (Evenden 1991). It is defined by 



Re = -2- 



180° /ln(cos£) m(cos<f b ) 



tan£ 



tan 2 £ b 



tan<f 



(71) 



where 

ft = 



90° -8 

2 ' 
90° - Bb 



When 8b approaches 90°, the second term of Eq. (71) ap- 
proaches its asymptotic value of -i. For all 8b, the projection 
is unbounded at the native south pole. Inversion of Eq. (71), a 
transcendental equation in 8, must be done via iterative meth- 
ods. 

The FITS keyword PVLla, attached to latitude coordi- 
nate i, will be used to specify 8b in degrees with a default of 
90°. This projection is illustrated in Fig. 14. 

5.2. Cylindrical projections 

Cylindrical projections are so named because the surface of 
projection is a cylinder. The native coordinate system is chosen 
to have its polar axis coincident with the axis of the cylinder. 
Meridians and parallels are mapped onto a rectangular gratic- 
ule so that cylindrical projections are described by formula? 
which return x and y directly. Since all cylindrical projections 
are constructed with the native coordinate system origin at the 
reference point, we set 



(00, #o)cylindrical - (0, 0) . 



(72) 



11 Devised in 1861 by astronomer royal George Biddell Airy, 1801- 
1892. 
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Fig. 15. Geometry of cylindrical projections. 
Furthermore, all cylindrical projections have 



(73) 



Cylindrical projections are often chosen to map the regions ad- 
jacent to a great circle, usually the equator, with minimal dis- 
tortion. 

5.2.1. CYP: cylindrical perspective 

Figure 15 illustrates the geometry for the construction of cylin- 
drical perspective projections. The sphere is projected onto a 
cylinder of radius A spherical radii from points in the equato- 
rial plane of the native system at a distance fi spherical radii 
measured from the center of the sphere in the direction oppo- 
site the projected surface. The cylinder intersects the sphere at 
latitudes 6 X = cos -1 A. It is straightforward to show that 



x = A0, 

y 



180°/ fi + A \ . 

= — sinfl. 

n \yU + cos#/ 



This may be inverted as 

x 



8 = arg(l,77) + sin 



where 



n 



y 



180° n + A 



(74) 
(75) 

(76) 
(77) 

(78) 



Note that all values of \x are allowable except fi = -A. For 
FITS purposes, we define the keywords PVLla to convey y. 
and PV/_2a for A, both measured in spherical radii, both with 
default value 1, and both attached to latitude coordinate i. 

The case with /j. = oo is covered by the class of cylindrical 
equal area projections. No other special-cases need be defined 
since cylindrical perspective projections have not previously 
been used in FITS. Aliases for a number of special cases are 
listed in Appendix C, Table A. 1 . Probably the most important 



Fig. 16. Gall's stereographic projection, CYP with yL = 1, # x = 45°; no 
limits. 
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Fig. 17. Lambert's equal area projection, CEA with A = 1; no limits. 

of these is Gall's stereographic projection, which minimizes 
distortions in the equatorial regions. It has jj. = I, A = V2/2, 
giving 

x = 0/V2, 

180° 



y = 



1 + 



It is illustrated in Fig. 16. 



5.2.2. CEA: cylindrical equal area 

The cylindrical equal area projection is the special case of the 
cylindrical perspective projection with // = 00. It is conformal 
at latitudes +9 C where A = cos 2 6 C . The formulae are 



x = (/>, 

180° sin 6 

which reverse to 



q> = x, 
9 = sin~ 



\180° / 



(79) 
(80) 

(81) 
(82) 



Note that the scaling parameter A is applied to the y-coordinate 
rather than the x-coordinate as in CYP. The keyword PVLla 
attached to latitude coordinate i is used to specify the dimen- 
sionless A with default value 1 . 

Lambert's 12 equal area projection, the case with A = 1, is 
illustrated in Fig. 17. It shows the extreme compression of the 

12 The mathematician, astronomer and physicist Johann Heinrich 
Lambert (1728-1777) was the first to make significant use of calcu- 
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Fig. 18. The plate carree projection (CAR); no limits. 

parallels of latitude at the poles typical of all cylindrical equal 
area projections. 
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5.2.3. CAR: plate carree 



Fig. 19. Mercator's projection (MER); diverges at 9 = ±90°. 



The equator and all meridians are correctly scaled in the plate 

carree projection 13 , whose main virtue is that of simplicity. Its This projection, illustrated in Fig. 19, has been widely used 
formulae are m navigation since it has the property that lines of constant 

bearing (known as rhumb lines or loxodromes) are projected as 
x = <f> , (83) straight lines. This is a direct result of its conformality and the 

y — q (34) fact that its meridians do not converge. 

Refer to Sect. 6.1.4 for a discussion of the usage of MER in 
The projection is illustrated in Fig. 18. AIPS. 



5.2.4. MER: Mercator 



5.3. Pseudocylindrical and related projections 



Since the meridians and parallels of all cylindrical projections 
intersect at right angles the requirement for conformality re- 
duces to that of equiscaling at each point. This is expressed by 
the differential equation 



dy - 1 dx 
89 ~ cos 0d~4> 7 

the solution 14 of which gives us Mercator's projection: 



<f>, 
180° 



y = 



■ In tan 



90° + 6 



with inverse 



(9 = 2 tan" 1 (e^ /180 °) - 90° . 



(85) 

(86) 
(87) 

(88) 



lus in constructing map projections. He formulated and gave his name 
to a number of important projections as listed in Table A.l. 

13 Although colloquially referred to as "Cartesian", 
Claudius Ptolemy (ca. 90-ca. 170 a.d.), influential cartographer 
and author of the Ptolemaic model of the solar system, credits 
Marinus of Tyre with its invention in about a.d. 100, thus predating 
Descartes by some 1500 years. 

14 Gerardus Mercator (1512-1594), a prominent Flemish map- 
maker, effectively solved this equation by numerical integration. 
Presented in 1569 it thus predates Newton's theory of fluxions by 
nearly a century. 



Pseudocylindricals are like cylindrical projections except 
that the parallels of latitude are projected at diminishing 
lengths towards the polar regions in order to reduce lat- 
eral distortion there. Consequently the meridians are curved. 
Pseudocylindrical projections lend themselves to the con- 
struction of interrupted projections in terrestrial cartography. 
However, this technique is unlikely to be of use in celestial 
mapping and is not considered here. Like ordinary cylindri- 
cal projections, the pseudocylindricals are constructed with 
the native coordinate system origin at the reference point. 
Accordingly we set 



(<f>0, 0o)pseudocylindrical - (0, 0) . 



(89) 



The Hammer-Aitoff projection is a modified zenithal projec- 
tion, not a pseudocylindrical, but is presented with this group 
on account of its superficial resemblance to them. 

5.3.1. SFL: Sanson-Flamsteed 

Bonne's projection (Sect. 5.5.1) reduces to the pseudocylindri- 
cal Sanson-Flamsteed 15 projection when 0\ = 0. Parallels are 
equispaced and projected at their true length which makes it an 



15 Nicolas Sanson d' Abbeville (1600-1667) of France and John 
Flamsteed (1646-1719), the first astronomer royal of England, pop- 
ularized this projection, which was in existence at least as early as 
1570. 
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Fig. 20. Sanson-Flamsteed projection (SFL); no limits. 




Fig. 21. Parabolic projection (PAR); no limits. 



equal area projection. The formulae are 
x = 0cos#, 

y = e, 

which reverse into 



cosy 

y- 



(90) 
(91) 

(92) 
(93) 



This projection is illustrated in Fig. 20. Refer to Sect. 6.1.4 for 
a discussion relating SFL to the GLS projection in AIPS. 

5.3.2. PAR: parabolic 

The parabolic or Craster pseudocylindrical projection is illus- 
trated in Fig. 21. The meridians are projected as parabolic arcs 
which intersect the poles and correctly divide the equator, and 
the parallels of latitude are spaced so as to make it an equal 
area projection. The formulae are 



<*|2cosy-l| . 



y = 180° sin 



6 



with inverse 
6 = 3 sin 
180 c 



1 180°/ 



<P = 



n l-4(y/180°) 



• 



(94) 
(95) 

(96) 
(97) 




Fig. 22. Mollweide's projection (M0L); no limits. 




Fig. 23. Hammer- Aitoff projection (AIT); no limits. 



5.3.3. M0L: Mollweide's 

In Mollweide's pseudocylindrical projection 16 , the meridians 
are projected as ellipses that correctly divide the equator and 
the parallels are spaced so as to make the projection equal area. 
The formulae are 



x = 0cosy, 

n 

r-180° . 

y = V2 sin y, 



(98) 
(99) 



where y is defined as the solution of the transcendental equation 

(100) 



7 sin 27 
sine = + - 



90° n 

The inverse equations may be written directly without refer 
ence to 7, 



7TX/\2J2-( 



: y)2 



180' 



6 = sin 



1 

90° 

y 



n y 



180 



180° V2V 

°a/ 2_ W )2 )' 



(101) 



(102) 



Mollweide's projection is illustrated in Fig. 22. 

16 Presented in 1805 by astronomer and mathematician Karl 
Brandan Mollweide (1774-1825). 
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Fig. 24. Construction of conic perspective projections (COP) and the resulting graticules; (left) two-standard projection with 0\ = 20°, 82 = 10°; 
(right) one-standard projection with 61 = d 2 = 45°. Both projections have Q d = 45° and this accounts for their similarity. Both diverge at 
9 = 6» a ± 90°. 



5.3.4. AIT: Hammer-Aitoff 



where 



The Hammer-Aitoff 17 projection illustrated in Fig. 23 is devel- 
oped from the equatorial case of the zenithal equal area projec- 
tion by doubling the equatorial scale and longitude coverage. 
The whole sphere is mapped thereby while preserving the equal 
area property. Note, however, that the equator is not evenly di- 
vided. 

This projection reduces distortion in the polar regions com- 
pared to pseudocylindricals by making the meridians and paral- 
lels more nearly orthogonal. Together with its equal area prop- 
erty this makes it one of most commonly used all-sky projec- 
tions. 

The formulae for the projection and its inverse are derived 
in Greisen (1986) and Calabretta (1992) among others. They 
are 



180° 



x = 2ycos0sin — , 
y = ysinfl, 



(103) 
(104) 



17 David Aitoff (1854-1933) developed his projection from the 
zenithal equidistant projection in 1889 and in 1892 Ernst Hammer 
(1858-1925) applied his idea more usefully to the zenithal equal area 
projection. See Jones (1993). 



n V 1 + cos6>cos(0/2) ' 
The reverse equations are 

^ 2arg ( 2z2 - 1 'T^r)' 

^sm'(^Z), 
where 

z= Ji-(—-) 2 -(— y -)\ 

V 1 180° 4/ 1 180° 2/ 
= +cos6»cos^j. 

that i 

the usage of AIT in MPS. 



(105) 

(106) 
(107) 

(108) 
(109) 



Note that \ <Z 2 < \. Refer to Sect. 6.1.4 for a discussion of 



5.4. Conic projections 

In conic projections the sphere is thought to be projected onto 
the surface of a cone which is then opened out. The native co- 
ordinate system is chosen so that the poles are coincident with 
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the axis of the cone. Native meridians are then projected as uni- 
formly spaced rays that intersect at a point (either directly or by 
extrapolation), and parallels are projected as equiangular arcs 
of concentric circles. 

Two-standard conic projections are characterized by two 
latitudes, 8\ and 82, whose parallels are projected at their true 
length. In the conic perspective projection these are the lati- 
tudes at which the cone intersects the sphere. One-standard 
conic projections have 0\ = 62 and the cone is tangent to the 
sphere as shown in Fig. 24. Since conies are designed to mini- 
mize distortion in the regions between the two standard paral- 
lels they are constructed so that the point on the prime meridian 
mid-way between the two standard parallels maps to the refer- 
ence point so we set 



(00, #0)conic = (0, 6> a ) , 

where 

a = (0i+0 2 )/2. 



(HO) 



(111) 



Being concentric, the parallels may be described by Rg, the ra- 
dius for latitude 8, and A , the angle for longitude (p. An offset 
in y is also required to force (x,y) = (0, 0) at ((/>, 6) = (0, a ). 
All one- and two-standard conies have 



C<f>, 



(112) 



where C, a constant known as the constant of the cone, is such 
that the apical angle of the projected cone is 360° C. Since the 
standard parallels are projected as concentric arcs at their true 
length we have 



_ 180° cos 0! _ 180° cos 2 
nRg l nRg 2 

Cartesian coordinates in the plane of projection are 

x - R g sin(C(f)), 

y = -R e cos(C<f>) + Y , 

and these may be inverted as 
Re = sign0 aA /x 2 + (y o -)O 2 , 



(113) 



(114) 
(115) 



(116) 
(117) 



To complete the inversion the equation for 8 as a function of R e 
is given for each projection. The equations given here correctly 
invert southern conies, i.e. those with d < 0. An example is 
shown in Fig. 25. 

The conies will be parameterized in FITS by a (given by 
Eq. (Ill)) and 77 where 



17 = |0! - 6 2 \/2 . 



(118) 



The keywords PVLla and PV;'_2a attached to latitude coordi- 
nate i will be used to specify a and 77 respectively, both mea- 
sured in degrees. PVLla has no default while PV/_2a defaults 
to 0. It is recommended that both keywords always be given. 




Fig. 25. Conic equal area projection (COE) with 61 = -20°, and 2 = 
-70°, also illustrating the inversion of southern hemisphere conies; no 
limits. 

Where 6\ and 02 are required in the equations below they may 
be computed via 

9i = 8,-11, (H9) 
82 = e a + n- (120) 

This sets 82 to the larger of the two. The order, however, is 
unimportant. 

As noted in Sect. 5, the zenithal projections are special 
cases of the conies with 8\ — 82 = 90°. Likewise, the cylin- 
drical projections are conies with 6\ = -82- However, we 
strongly advise against using conies in these cases for the rea- 
sons given previously. Nevertheless, the only formal require- 
ment on 8\ and 82 in the equations presented below is that 
-90° < 01,02 < 90°. 



5.4.1 . COP: conic perspective 

Development of Colles' conic perspective projection is shown 
in Fig. 24. The projection formulae are 



sin 8. d , 
180° 



• cos rj [cot a - tan(0 - a )] , 



n 
180° 

Y = cos 77 cot 8. d . 

n 



The inverse may be computed with 

Jj_ Re_ 

180° cost; 



' = a + tan 1 |cot0 a - 



(121) 
(122) 

(123) 
(124) 



5.4.2. COE: conic equal area 

The standard parallels in Alber's conic equal area projection are 
projected as concentric arcs at their true length and separated 
so that the area between them is the same as the corresponding 
area on the sphere. The other parallels are then drawn as con- 
centric arcs spaced so as to preserve the area. The projection 
formula? are 



y/2, 



(125) 



20 
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Fig. 27. Conic orthomorphic projection (COO) with 6\ = 20° and 6 2 
70°; diverges ate = -90°. 



Fig. 26. Conic equidistant projection (COD) with 6\ = 20° and 8 2 
70° ; no limits. 



180° 2 

Rg = yfl + sin 8\ sin 62 - y sin , 



Ya = 



n y 
180° 2 



n y 



yjl +sin6»i sin6» 2 -ysin((6»i +<9 2 )/2), (127) 



where 

y = sin 0\ + sin 62 . 

The inverse may be computed with 



(128) 



= sin 



(129) 



_, / 1 sin 0i sin 02 / 

in - + y 

\y y y \360° 

This projection is illustrated in Fig. 25. 
5.4.3. COD: conic equidistant 

In the conic equidistant projection the standard parallels are 
projected at their true length and at their true separation. The 
other parallels are then drawn as concentric arcs spaced at their 
true distance from the standard parallels. The projection for- 
mulae are 

180° sin a sin 77 



n Tj 

Rg = #a - 9 + J] COt T] cot # a , 
Yq = T] COt 77 COt 6 a . 

The inverse may be computed with 
= 9-d + i] cot T] cot G a - Rg. 
For 6\ = 02 these expressions reduce to 
C = sin 6. d , 

180° 

Rg = 0^-0+ COt# a , 

n 

180° 

Yq = cot0 a , 



(130) 

(131) 
(132) 

(133) 

(134) 
(135) 

(136) 



and inverse 

180° 

= a + cot 0- d -Rg. 



(137) 



(126) -p n j s p ro j ec ti on [ s illustrated in Fig. 26. 



5.4.4. COO: conic orthomorphic 

The requirement for conformality of conic projections is 

8Rg = -nRg 
80 



-C. 



(138) 



180° cos0 

Solution of this differential equation gives rise to the formula; 
for Lambert's conic orthomorphic projection: 



C 



V cos 9i I 



Rg = 

Yo = 
where 



In 

<!> 

<!> 

180° 



tan 


; 9o vi 


tan 


^9O°-0 a j 





cos# 



C 



180° 



tan (22^) 

COS 62 



tan 



(139) 

(140) 
(141) 

(142) 
(143) 



The inverse may be computed with 



= 90° -2 tan" 1 



-A 



(144) 



When 0\ = 82 the expression for C may be replaced with C = 
sin 9\. This projection is illustrated in Fig. 27. 
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Fig. 28. Bonne's projection (BON) with di = 45°; no limits. 




5.5. Polyconic and pseudoconic projections 

Polyconics are generalizations of the standard conic projec- 
tions; the parallels of latitude are projected as circular arcs 
which may or may not be concentric, and meridians are curved 
rather than straight as in the standard conies. Pseudoconics are 
a sub-class with concentric parallels. The two polyconics pre- 
sented here have parallels projected at their true length and use 
the fact that for a cone tangent to the sphere at latitude 6\, as 
shown in Fig. 24, we have R 6l = cot 6\ . Since both are 
constructed with the native coordinate system origin at the ref- 
erence point we set 



(00, #o)polyconic = (0, 0) . 



5.5.1 . BON: Bonne's equal area 



(145) 



In Bonne's pseudoconic projection 18 all parallels are projected 
as concentric equidistant arcs of circles of true length and true 
spacing. This is sufficient to guarantee that it is an equal area 
projection. It is parameterized by the latitude 6\ for which 
R 8l = i^Q- cot 6 l . The projection is conformal at this latitude 
and along the central meridian. The equations for Bonne's pro- 
jection become divergent for 8\ = and this special case is 
handled as the Sanson-Flamsteed projection. The projection 
formulae are 



y = -R g cos A + Y , 



where 



180° 



* nRg 
Re = Y Q - 
Y - 18 ° C 

7T 



0COS0, 



COt Bi+6 



(146) 
(147) 



(148) 
(149) 
(150) 



18 The ancestry of this important projection may be traced back to 
Ptolemy. Attribution for its invention is uncertain but it was certainly 
used before the birth of Rigobert Bonne (1727-1795) who did much 
to popularize it. 



Fig. 29. Polyconic projection (PCO); no limits. 



The inverse formula? are then 



= Yq - Rg , 

where 



Rg = sign8 l ^x 2 + (Y Q -y) 2 , 
'Y () -y x 



^ = arg 



Rg Rf, 



(151) 
(152) 



(153) 
(154) 



This projection is illustrated in Fig. 28. The keyword PVLla at- 
tached to latitude coordinate i will be used to give 8\ in degrees 
with no default value. 



5.5.2. PCO: polyconic 

Each parallel in Hassler's polyconic projection is projected 
as an arc of a circle of radius cot 6 at its true length, 
360° cos 8, and correctly divided. The scale along the central 
meridian is true and consequently the parallels are not concen- 
tric. The projection formula? are 



180° 

x = cot 8 sm((f> sin 8) , 

n 

180° 

y = 8 + cot 8 [1 - cos(0 sin 8)] . 

n 

Inversion requires iterative solution for 8 of the equation 

, 360° , 
x 2 (y - 8) cot6» + {y-df =0. 

7T 

Once 8 is known is given by 

1 /180° \ 

d> = arg (y - 8) tan 8, x tan 8\ . 

sm8 \ n ) 

The polyconic projection is illustrated in Fig. 29. 



(155) 
(156) 

(157) 

(158) 
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5.6. Quad-cube projections 

Quadrilateralized spherical cube (quad-cube) projections be- 
long to the class of polyhedral projections 19 in which the sphere 
is projected onto the surface of an enclosing polyhedron, typi- 
cally one of the five Platonic polyhedra: the tetrahedron, hexa- 
hedron (cube), octahedron, dodecahedron, and icosahedron. 

The starting point for polyhedral projections is the 
gnomonic projection of the sphere onto the faces of an enclos- 
ing polyhedron. This may then be modified to make the pro- 
jection equal area or impart other special properties. However, 
minimal distortion claims often made for polyhedral projec- 
tions should be balanced against the many interruptions of the 
flattened polyhedron; the breaks should be counted as extreme 
distortions. Their importance in astronomy is not so much in 
visual representation but in solving the problem of distributing 
N points as uniformly as possible over the sphere. This may 
be of particular importance in optimizing computationally in- 
tensive applications. The general problem may be tackled by 
pixelizations such as HEALPix (Hierarchical Equal Area iso- 
Latitude PIXelization, Gorski 1999), which define a distribu- 
tion of points on the sphere but do not relate these to points in 
a plane of projection and are therefore outside the scope of this 
paper. 

The quad-cubes have been used extensively in the COBE 
project and are described by Chan & O'Neill (1975) and 
O'Neill & Laubscher (1976). The icosahedral case has also 
been studied by Tegmark (1996). It is close to optimal, provid- 
ing a 10% improvement over the cubic case. However, we have 
not included it here since it relies on image pixels being or- 
ganized in an hexagonal close-packed arrangement rather than 
the simple rectangular arrangement supported by FITS. 

The six faces of quad-cube projections are numbered and 
laid out as 



4 3 2 1 4 3 2 
5 

where faces 2, 3 and 4 may appear on one side or the other 
(or both). The layout used depends only on the FITS writer's 
choice of C in Table 3. It is also permissible to split faces be- 
tween sides, for example to put half of face 3 to the left and 
half to the right to create a symmetric layout. FITS readers 
should have no difficulty determining the layout since the ori- 
gin of (x, y) coordinates is at the center of face 1 and each face 
is 90° on a side. The range of x therefore determines the lay- 
out. Other arrangements are possible and Snyder (1993) illus- 
trates the "saw-tooth" layout of Reichard's 1803 map of the 
Earth. While these could conceivably have benefit in celestial 
mapping we judged the additional complication of representing 
them in FITS to be unwarranted. The layout used in the COBE 
project itself has faces 2, 3, and 4 to the left. 

19 Polyhedral projections date from renaissance times when the artist 
and mathematician Albrecht Diirer (1471-1528) described, although 
did not implement, the tetrahedral, dodecahedral, and icosahedral 
cases. Snyder (1993) traces subsequent development into the twentieth 
century, including one by R. Buckminster Fuller onto the non-Platonic 
"cuboctahedron" with constant scale along each edge. 



Table 3. Assignment of parametric variables and central longitude and 
latitude by face number for quadrilateralized spherical cube projec- 
tions. 



Face f 7] £ (p c C 






m 


-I 


n 




0° 




90° 


1 


m 


n 


I 




0° 




0° 


2 


-I 


n 


m 


- 270° 


or 


90° 


0° 


3 


-m 


n 


-I 


-180° 


or 


180° 


0° 


4 


I 


n 


-m 


-90° 


or 


270° 


0° 


5 


in 


I 


-n 




0° 




-90° 



The native coordinate system has its pole at the center of 
face and origin at the center of face 1 (see Fig. 30) which is 
the reference point whence 

(0O,#o) q u a d-cube = (O,O). (159) 

The face number may be determined from the native longitude 
and latitude by computing the direction cosines: 

I = COS0COS0, 

m = cos#sin0, (160) 
n = sin#. 

The face number is that which maximizes the value of £ in 
Table 3. That is, if f is the largest of n, I, m, -I, -m, and -n, 
then the face number is through 5, respectively. Each face 
may then be given a Cartesian coordinate system, (£, rj), with 
origin in the center of each face as per Table 3. The formulae 
for quad-cubes are often couched in terms of the variables 

X = (161) 
<A = tt/(. (162) 

Being composed of six square faces the quad-cubes admit 
the possibility of efficient data storage in FITS. By stacking 
them on the third axis of a three-dimensional data structure the 
storage required for an all-sky map may be halved. This axis 
will be denoted by a CTYPE/a value of 'CUBEFACE'. In this case 
the value of (x,y) computed via Eq. (1) for the center of each 
face must be (0, 0) and the FITS interpreter must increment this 
by ((f> c ,6 c ) using its choice of layout. Since the CUBEFACE axis 
type is purely a storage mechanism the linear transformation of 
Eq. (1) must preserve the CUBEFACE axis pixel coordinates. 

5.6.1 . TSC: tangential spherical cube 

While perspective quad-cube projections could be developed 
by projecting a sphere onto an enclosing cube from any point 
of projection, inside or outside the sphere, it is clear that only 
by projecting from the center of the sphere will every face be 
treated equally. Thus the tangential spherical cube projection 
(TSC) consists of six faces each of which is a gnomonic pro- 
jection of a portion of the sphere. As discussed in Sect. 5.1.1, 
gnomonic projections map great circles as straight lines but un- 
fortunately diverge very rapidly away from the poles and can 
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Fig. 30. Tangential spherical cube projection (TSC); no limits. 

only represent a portion of the sphere without extreme distor- 
tion. The TSC projection partly alleviates this by projecting 
great circles as piecewise straight lines. To compute the for- 
ward projection first determine x an d ^ as described above, 
then 



x = 4> C +45° X , 
y = 6» c +45>. 



(163) 
(164) 



To invert these first determine to which face the (x, y) coordi- 
nates refer, then compute 



X = (*-0 c )/45°, 
<A = (y-e c )/45°, 

then 

i = n^\+x 2 +* 2 . 

Once £ is known £ and r\ are obtained via 

£ = *f. 

T] = 



(165) 
(166) 



(167) 



(168) 
(169) 



The direction cosines (/, m, n) may be identified with (£, jj, £ ) 
with with the aid of Table 3, whence (0, ff) may readily be com- 
puted. The projection is illustrated in Fig. 30 for the full sphere. 

5.6.2. CSC: COBE quadrilateralized spherical cube 

The COBE quadrilateralized spherical cube projection illus- 
trated in Fig. 31 modifies the tangential spherical cube projec- 
tion in such a way as to make it approximately equal area. The 
forward equations are 



x = <p c +45°F(x,i(f), 
y = c +45°F(ifr,x), 

where the function F is given by 



(170) 
(171) 



Fig. 31. COBE quadrilateralized spherical cube projection (CSC); no 
limits. 



^ 2 d 



-x 2 ) 



T + (M - T)x 2 + 



i=0 7=0 



* 3 (i -x 2 ) 



(172) 



dj and Di are derived from c*. and d* as given by Chan & 
O'Neill (1975). The other parameters are given by exact for- 
mula? developed by O'Neill & Laubscher (1976), who provide 
the numeric values of their parameters in tables and software 
listings. Both disagree with their formulae, but the software list- 
ings do contain the actual numeric parameters still in use for the 
COBE Project (Immanuel Freedman, private communication, 
1993). They are 



7* = 


1.37484847732 


M = 


0.004869491981 


r = 


-0.13161671474 


Qi = 


-0.159596235474 


Coo = 


0.141189631152 


Cio = 


0.0809701286525 


Coi = 


-0.281528535557 


C20 = 


-0.178251207466 


C n = 


0.15384112876 


Cq2 = 


0.106959469314 


Do = 


0.0759196200467 


Di = 


-0.0217762490699 



Chan & O'Neill (1975) actually defined the projection via 
the inverse equations. Their formulation may be rewritten in a 
more convenient form which is now the current usage in the 
COBE Project (Immanuel Freedman, private communication, 
1993): 



X = f(x-<f> c ,y-e c ), 
4* = f(y-0c,x-4> c ), 



(173) 
(174) 
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Fig. 32. Quadrilateralized spherical cube projection (QSC); no limits. 



where 



f(x-<f> c ,y-0 c ) = 

N N-j 

X + X ( 1 - X 2 ) P u X 2i Y 2j , 

7=0 i=0 

and 



(175) 



X = (x-<t> c )/45°, 
Y = (y- 9 C )IA5° . 



For COBE, N 
to be 



6 and the best-fit parameters have been taken 



^00 




-0.27292696 


Po4 




0.93412077 


Pw 




-0.07629969 


P50 




0.25795794 


Poi 




-0.02819452 


P41 




1.71547508 


P20 




-0.22797056 


P32 




0.98938102 


Pn 




-0.01471565 


P23 




-0.93678576 


P02 




0.27058160 


Pu 




-1.41601920 


P30 




0.54852384 


P05 




-0.63915306 


Pll 




0.48051509 


P&Q 




0.02584375 


Pn 




-0.56800938 


PSI 




-0.53022337 


P03 




-0.60441560 


P42 




-0.83180469 


P40 




-0.62930065 


P33 




0.08693841 


P31 




-1.74114454 


P24 




0.33887446 


Pl2 




0.30803317 


Pl5 




0.52032238 


Pn 




1.50880086 


P06 




0.14381585 



Given the face number, x, and if/, the native coordinates ((/>, 6) 
may be computed as for the tangential spherical cube projec- 
tion. 

Equations (172) and (175), the forward and reverse projec- 
tion equations used by COBE, are not exact inverses. Each set 
could of course be inverted to any required degree of precision 
via iterative methods (in that case Eq. (175) should be taken to 
define the projection). However, the aim here is to describe the 
projection in use within the COBE project. One may evaluate 
the closure error in transforming (x - <f> c ,y - 6 C ) to iff) with 



Eq. (175) and then transforming back to (x 
Eq. (172), i.e. 



0c, y - 0c) with 



Ei 



(FWwjlfiy^XiV-x,) 2 
+ (F(f(y j ,x i ),f(x i ,y j ))-y j f . 



The COBE parameterization produces an average error of 4.7 
arcsec over the full field. The root mean square and peak errors 
are 6.6 and 24 arcsec, respectively. In the central parts of the 
image \Y\ < 0.8), the average and root mean square errors 
are 5.9 and 8.0 arcsec, larger than for the full field. 

Measures of equal-area conformance obtained for Eq. (175) 
show that the rms deviation is 1.06% over the full face and 
0.6% over the inner 64% of the area of each face. The maxi- 
mum deviation is +13.7% and -4.1% at the edges of the face 
and only ±1 .3% within the inner 64% of the face. 

5.6.3. QSC: quadrilateralized spherical cube 

O'Neill & Laubscher (1976) derived an exact expression for 
an equal-area transformation from a sphere to the six faces 
of a cube. At that time, their formulation was thought to be 
computationally intractable, but today, with modern computers 
and telescopes of higher angular resolution than COBE, their 
formulation has come into use. Fred Patt (1993, private com- 
munication) has provided us with the inverse of the O'Neill 
& Laubscher formula and their expression in Cartesian coordi- 
nates. 

O'Neill & Laubscher's derivation applies only in the quad- 
rant -45° < < 45° and must be reflected into the other 
quadrants. This has the effect of making the projection non- 
differentiable along the diagonals as is evident in Fig. 32. To 
compute the forward projection first identify the face and find 
(£, 77, and (</> c , C ) from Table 3. Then 



(x,y) = (0c,0 c ) + 



where 



(11, v) if\{\>\n\ 

(v,u) otherwise 



u = 45°S 
u 

V= 15= 



1-f 



\ 1 - 1/ y/2 + co 

tan -1 (oj) - sin -1 



W2(l+^ 2 ) 



(176) 

(177) 
(178) 



a> = 



I if£|>|77l 

1 otherwise ' 

+ 1 iff > \tj\ oxt 1 > |f| 

-1 otherwise 



To compute the inverse first identify the face from the (x,y) 
coordinates, then determine (u, v) via 



(u,v) 



Then 



(x -(p c ,y- 0c) if \x - <p c \ >\y-0 c 
(y — 6 c ,x- C ) otherwise 



^ 1 -(^) 2 



1 - 



V2 



(179) 



(180) 
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where 



sin(15°v/«) 



cos(15°v/m)- 1/ V2 
If |x-0 c | > |y-0 c |then 



otherwise 



V 1 + to 2 ' 



n = 



1 +C0 2 ' 



(181) 



(182) 
(183) 

(184) 
(185) 



Given the face number and (£, 77, £), the native coordinates (0, 0) 
may be computed with reference to Table 3 as for the tangential 
spherical cube projection. 

6. Support for the AIPS convention 

A large number of FITS images have been written using the 
AIPS coordinate convention and a substantial body of software 
exists to interpret it. Consequently, the AIPS convention has 
acquired the status of a de facto standard and FITS interpreters 
will need to support it indefinitely in order to obey the maxim 
"once FITS always FITS". Translations between the old and 
new system are therefore required. 



6.1. Interpreting old headers 

In the AIPS convention, CROTA; assigned to the latitude axis 
was used to define a bulk rotation of the image plane. Since this 
rotation was applied after CDELT; the translation to the current 
formalism follows from 



CDELT 1 
CDELT2 



(186) 



PCl.l PC1_2 
PC2_1 PC2_2 

cosp - sinp \ / CDELT1 
sinp cospjy CDELT2 

where we have used subscript 1 for the longitude axis, 2 for 
latitude, and written p for the value of CR0TA2. Equation (186), 
which includes the added constraint of preserving CDELT/ in 
the translation, is readily solved for the elements of the PC Lja 
matrix 



/ pci_: 
\ PC2_: 

where 



PC1_2 
PC2_2 



cosp -A sinp' 
j sin p cos p 



A = CDELT2/CDELT1 . 



(187) 



(188) 



Note that Eq. (187) defines a rotation if and only if A = ±1, 
which is often the case. In fact, the operations of scaling and 
rotation are commutative if and only if the scaling is isotropic, 
i.e. A = +1; for A = -1 the direction of the rotation is re- 
versed. However, whatever the value of A, the interpretation of 
Eq. (186) as that of a scale followed by a rotation is preserved. 



The translation for CDi_ja is simpler, effectively because 
the CDELT / have an implied value of unity and the constraint 
on preserving them in the translation is dropped: 



CD 1.1 CD1_2 
CD2_1 CD2_2 



CDELT 1 cosp -CDELT2 sinp' 
CDELT 1 sinp CDELT2 cosp 



(189) 



The expressions in Hanisch & Wells (1988) and Geldzahler 
& Schlesinger (1997) yield the same results as Eq. (189) for the 
usual left-handed sky coordinates and right-handed pixel coor- 
dinates, but can lead to an incorrect interpretation (namely, pos- 
sible sign errors for the off-diagonal elements) for other config- 
urations of the coordinate systems. The Hanisch & Wells draft 
and the coordinate portions of Geldzahler & Schlesinger are 
superseded by this paper. 

6.1.1. SIN 

The SIN projection defined by Greisen (1983) is here general- 
ized by the addition of projection parameters. However, these 
parameters assume default values which reduce to the simple 
orthographic projection of the AIPS convention. Therefore no 
translation is required. 

6.1.2. NCP 

The "north-celestial-pole" projection defined by Greisen 
(1983) is a special case of the new generalized SIN projection. 
The old header cards 

CTYPE1 = 'RA NCP' , 

CTYPE2 = 'DEC- -NCP' , 

should be translated to the current formalism as 

CTYPE1 = 'RA---SIN' , 

CTYPE2 = 'DEC--SIN' , 

PV2_1 = 0, 

PV2_2 = cot5 . 

6.1.3. TAN, ARC and STG 

The TAN, ARC and STG projections defined by Greisen (1983, 
1986) are directly equivalent to those defined here and no trans- 
lation is required. 

6.1.4. AIT, GLS and MER 

Special care is required in interpreting the AIT (Hammer- 
Aitoff), GLS (Sanson-Flamsteed), and MER (Mercator) projec- 
tions in the AIPS convention as defined by Greisen (1986). As 
explained in Sect. 7.1, the AIPS convention cannot represent 
oblique celestial coordinate graticules such as the one shown 
in Fig. 2. CRVAL/ for these projections in AIPS does not corre- 
spond to the celestial coordinates (ao, 6q) of the reference point, 
as understood in this formalism, unless they are both zero in 
which case no translation is required. 



26 



Mark R. Calabretta and Eric W. Greisen: Representations of celestial coordinates in FITS 



A translation into the new formalism exists for non-zero 
CRVAL ;' but only if CROTA; is zero. It consists of setting CRVAL ;' 
to zero and adjusting CRPIX j and CDELT; accordingly in the 
AIPS header whereupon the above situation is obtained. The 
corrections to CRPIX j are obtained by computing the pixel co- 
ordinates of (a, S) = (0, 0) within the AIPS convention. For AIT 
and MER (but not GLS), CDELT; must also be corrected for the 
scaling factors f a and f s incorporated into the AIPS projection 
equations. 

Of the three projections only GLS is known to have been 
used with non-zero CRVAL ;'. Consequently we have renamed it 
as SFL as a warning that translation is required. 

6.2. Supporting old interpreters 

As mentioned in Sect. 6, FITS interpreters will need to recog- 
nize the AIPS convention virtually forever. It stands to reason, 
therefore, that if modern FITS-writers wish to assist older FITS 
interpreters they may continue to write older style headers, as- 
suming of course that it is possible to express the coordinate 
system in the AIPS convention. 

Modern FITS-writers must not attempt to help older inter- 
preters by including CROTA; together with the new keyword 
values (assuming the combination of CDELT; and PCLja ma- 
trix, or CDLja matrix, is amenable to such translation). We 
make this requirement primarily to minimize confusion. 

Assuming that a header has been developed using the 
present formalism the following test may be applied to deter- 
mine whether the combination of CDELT ia and PC ;'_/'« matrix 
represents a scale followed by a rotation as in Eq. (186). Firstly 
write 

/CD 1.1 CD1_2 \ _ 
CD2_1 CD2_2 



CDELT 1 
CDELT2 



PCl.l PC1_2 
PC2_1 PC2_2 



(190) 



where 1 is the longitude coordinate and 2 the latitude coordi- 
nate, then evaluate p a and p b as 

farg( CD1_1, CD2_1) if CD2_1 > 
p. d = I if CD2.1 = , 

[ arg(-CDl_l,-CD2_l) ifCD2_l<0 

f arg(-CD2_2, CD1_2) if CD1_2 > 
p b = I if CD1.2 = . (191) 

[arg( CD2_2,-CD1_2) ifCDl_2<0 

If p a = pb to within reasonable precision (Geldzahler & 
Schlesinger, 1997), then compute 

P = (p a +Pb)/2 (192) 

as the best estimate of the rotation angle, the older keywords 
are 



CDELT1 = CDl_l/cosp, 
CDELT2 = CD2_2/cosp, 
CR0TA2 = p. 



(193) 



Note that the translated values of CDELT; in Eqs. (193) may 
differ from the starting values in Eq. (190). 



Solutions for CR0TA2 come in pairs separated by 180°. The 
above formulae give the solution which falls in the half-open 
interval [0, 180°). The other solution is obtained by subtracting 
180° from CR0TA2 and negating CDELT 1 and CDELT2. While 
each solution is equally valid, if one makes CDELT 1 < and 
CDELT2 > then it would normally be the one chosen. 

Of course, the projection must be one of those supported 
by the AIPS convention, which only recognizes SIN, NCP, TAN, 
ARC, STG, AIT, GLS and MER. Of these, we strongly recommend 
that the AIPS version of AIT, GLS, and MER not be written be- 
cause of the problems described in Sect. 6.1.4. It is interesting 
to note that a translation does exist for the slant orthographic 
(SIN) projection defined in Sect. 5.1.5 to the simple ortho- 
graphic projection of AIPS. However, we advise against such 
translation because of the likelihood of creating confusion and 
so we do not define it here. The exception is where the SIN 
projection may be translated as NCP as defined in Sect. 6.1.2. 

7. Discussion 

7.1. Oblique projections 

The term oblique projection is often used when a projection's 
native coordinate system does not coincide with the coordi- 
nate system of interest. Texts on terrestrial cartography often 
separately name and derive projection equations for particu- 
lar oblique projections. Thus the Cassini, Transverse Mercator, 
and Bartholomew nordic projections are nothing more than the 
plate carree, Mercator, and Mollweide projections in disguise. 

The view of obliqueness as being a property of a projection 
arose mainly because of the computational difficulty of produc- 
ing oblique graticules in the days before electronic computers. 
The particular aspects chosen were those for which geomet- 
rical construction was possible or for which the mathematical 
formulation had a simple form, and this tended to entrench par- 
ticular oblique projections as separate entities. In addition, ter- 
restrial longitude and latitude are so closely tied to the visual 
representation of the Earth's surface that it is shown predomi- 
nantly in the usual north-south orientation. The only other natu- 
ral terrestrial coordinate system is that defined by the magnetic 
pole but the difference between it and the geographic graticule 
is sufficiently small that it is usually treated as a local correc- 
tion to the magnetic bearing. The special view of obliqueness 
was probably also reinforced by traditional methods of map 
making using sextant and chronometer, which were based on 
geographic longitude and latitude. 

The situation in astronomical cartography is quite differ- 
ent. The celestial sphere has a variety of natural coordinate 
systems - equatorial, ecliptic, galactic, etc. - and oblique and 
non-oblique graticules are often plotted together on the same 
map. It wouldn't make sense to describe such a projection as 
being simultaneously oblique and non-oblique; clearly it is the 
particular coordinate graticule which may be oblique, not the 
projection. Visually, an area of the sky may be seen in differ- 
ent orientations depending on whether it's rising, transiting, or 
setting, or whether seen from the northern or southern hemi- 
sphere. Moreover, oblique graticules arise as a normal feature 
of observations in optical, infrared, radio, and other branches 
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of astronomy, just as they do now in terrestrial mapping based 
on aerial and satellite photography. The center of the field of 
view, wherever it may happen to be, typically corresponds to 
the native pole of a zenithal projection. Thus the aim of this 
paper has been to provide a formalism whereby obliquity may 
be handled in a general way. 

Figures 33 and 34 illustrate the same four oblique celes- 
tial graticules for the zenithal equal area, conic equidistant, 
Hammer-Aitoff, and COBE quadrilateralized spherical cube 
projections (at variable (x,y) scale). For the sake of intercom- 
parison, these graticules are defined in terms of the celestial 
coordinates of the native pole, (a p , 6 V ), together with p by set- 
ting (<f> , 0o) = (0, 90°) as described in Sect. 2.5. 

The first graticule, A, when compared to the non-oblique 
native coordinate graticules presented earlier for each pro- 
jection, illustrates the effect of changing 6 V (and hence <5o). 
Comparison of graticules A and B shows that changing a v 
(and hence ao) results in a simple change in origin of longi- 
tude. Graticules A, C, and D show the more interesting effect of 
varying p (conveyed by LONPOLEa). For the zenithal projec- 
tions the result is indistinguishable from a bulk rotation of the 
image plane, but this is not the case for any other class of pro- 
jection. This explains why the role of LONPOLEa was covered 
by CROTA; in the AIPS convention for the zenithal projections 
introduced by Greisen (1983), but why this does not work for 
any other class of projection. 

7.2. Choice of projection 

The projected meridians and parallels in Figs. 33 and 34 also 
serve to illustrate the distortions introduced by the various pro- 
jections. In particular, the quad-cube projection, while doing a 
good job within each face, is very distorted over the sphere as 
a whole, especially where the faces join, so much so that it is 
difficult even to trace the path of some of the meridians and 
parallels. However, as indicated previously, this projection is 
designed for efficient numerical computation rather than visual 
representation of the sphere. 

This leads us to the question of the choice of projection 
for a particular application. In some cases the projection is the 
natural result of the geometry of the observation and there is 
no choice. For example, maps produced by rotational synthe- 
sis radio telescopes are orthographic (SIN) projections with na- 
tive pole at the phase center of the observation. Photographic 
plates produced by Schmidt telescopes are best described by 
the zenithal equidistant (ARC) projection, while the field of view 
of other optical telescopes is closer to a gnomonic (TAN) pro- 
jection. On the other hand, the great circle scanning technique 
of the Sloan Digital Sky Survey produces a cylindrical pro- 
jection. Similarly a map of the surface of the Moon as it ap- 
pears from Earth requires a zenithal perspective (AZP) projec- 
tion with /i « -220. The same is true for spacecraft generated 
images of distant moons and planets. All-sky cameras should 
be well served by a ZPN projection with empirically determined 
parameters. 

Sometimes observational data will be regridded onto a pro- 
jection chosen for a particular purpose. Equal area projections 



are often used since they preserve surface density and allow 
integration, whether numerical or visual, to be performed with 
reasonable accuracy by summing pixel values. The Hammer- 
Aitoff (AIT) projection is probably the most widely used for 
all-sky maps. However, the Sanson-Flamsteed (SFL) may be 
preferred as being easier to take measurements from. The hum- 
ble plate carree (CAR) excels in this regard and may be consid- 
ered adequate, say, for mapping a few degrees on either side 
of the galactic plane. In general terms, zenithal projections are 
good for mapping the region in the vicinity of a point, often a 
pole; cylindrical projections are good for the neighborhood of a 
great circle, usually an equator; and the conies are suitable for 
small circles such as parallels of latitude. A conic projection 
might be a good choice for mapping a hemisphere; distortion 
at most points is reduced compared to a zenithal projection al- 
though at the expense of the break between native longitudes 
±180°. However, in some cases this break might be considered 
an extreme distortion which outweighs other reductions in dis- 
tortion and, depending on the application, may mandate the use 
of a zenithal projection. Cartographers typically favor conies 
and polyconics for "Australia-sized" regions of the sphere and 
at this they excel. Oblique forms of the zenithal equal-area pro- 
jection are also commonly used. Celestial applications might 
include large scale maps of the Magellanic clouds. Tables 3.1 
and 4.1 of Snyder (1993) summarize actual usage in a variety 
of 19th and 20th century atlases. 

As mentioned in Sect. 5, some projections are not scaled 
true at the reference point. In most cases this is deliberate, usu- 
ally because the projection is designed to minimize distortion 
over a wide area; the scale will be true at some other place 
on the projection, typically along a parallel of latitude. If the 
projection is required to have (a, 6) « (aso,So) + (x,y) then a 
number of projections will serve provided also that P takes its 
default value. The following are always scaled true at the refer- 
ence point: SZP, TAN, SIN, STG, ARC, ZEA, CAR, MER, BON, PCO, 
SFL, and AIT. The following are scaled true at the reference 
point for the particular conditions indicated: AZP (y = 0), ZPN 
(P = 0,Pi = 1), AIR (6 b = 90°), CYP (A = 1), CEA (A = 1), 
COP (77 = 0), COD (77 = 0), C0E (77 = 0), and COO (77 = 0). The 
following are never scaled true at the reference point: PAR (but 
is scaled true in x), MOL, TSC, CSC, and QSC; the latter three are 
equiscaled at the reference point. 

Of course CDELT/a may be used to control scaling of (x,y) 
with respect to (771,772)- Thus, for example, the plate carree pro- 
jection (CAR) also serves as the more general equirectangular 
projection. 

7.3. Header interpretation examples 

We now consider three examples chosen to illustrate how a 
FITS reader would interpret a celestial coordinate header. 

Example 1 is a simple header whose interpretation is quite 
straightforward. 

Example 2 is more complicated; it also serves to illustrate 
the WCS header cards for 1) image arrays in binary tables, and 
2) pixel lists. 



28 



Mark R. Calabretta and Eric W. Greisen: Representations of celestial coordinates in FITS 






B 






Fig. 34. Oblique (a, 6) graticules plotted for the Hammer- Aitoff (AIT) and COBE quadrilateralized spherical cube (CSC) projections using the 
same parameters as Fig. 33. 
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Table 4. FITS header for example 1. 



123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789 



NAXIS 


= 




4 


/ 


4-dimensional cube 


NAXIS1 


= 




512 


/ 


x axis (fastest) 


NAXIS2 


= 




512 


/ 


y axis (2nd fastest) 


NAXIS 3 


= 




196 


/ 


z axis (planes) 


NAXIS4 


= 




1 


/ 


Dummy to give a coordinate 


CRPIX1 


= 




256 


/ 


Pixel coordinate of reference point 


CDELT1 


= 




-0.083 


/ 


18.8 arcsec per pixel 


CTYPE1 


= 


'RA— TAN' 




/ 


Gnomonic projection 


CRVAL1 


= 




45.83 


/ 


RA at reference point 


CUNIT1 


= 


'deg 




/ 


Angles are degrees always 


CRPIX2 


= 




257 


/ 


Pixel coordinate of reference point 


CDELT2 


= 




0.003 


/ 


18.8 arcsec per pixel 


CTYPE2 


= 


'DEC— TAN' 




/ 


Gnomonic projection 


CRVAL2 


= 




63.57 


/ 


Dec at reference point 


CUNIT2 


= 


'deg 




/ 


Angles are degrees always 


CRPIX3 


= 




1 


/ 


Pixel coordinate of reference point 


CDELT3 


= 




7128.3 


/ 


Velocity increment 


CTYPE3 


= 


'VELOCITY' 




/ 


Each plane at a velocity 


CRVAL3 


= 




500000.0 


/ 


Velocity in m/s 


CUNIT3 




'm/s 




/ 


metres per second 








i 
i 


/ 
/ 


rlAcl LUUIUlILaLt: UI Iclcl cllLc pU-LILL 


CDELT4 






1 


/ 


Required here. 


CTYPE4 




' STOKES ' 




/ 


Polarization 


CRVAL4 






1 


/ 


Unpolarized 


CUNIT4 








/ 


Conventional unitless = I pol 


LONPOLE 






188 


/ 


Native longitude of celestial pole 


RADESYS 




' FK5 




/ 


Mean IAU 1984 equatorial coordinates 


EQUINOX 






2080.0 


/ 


Equator and equinox of J 2808. 8 



Example 3 highlights a subtle problem introduced into a 
header by the FITS writer and considers how this should be 
corrected. As such, it introduces the generally more difficult 
task of composing WCS headers, considered in greater detail 
in Sect. 7.4. 

7.3.1. Header interpretation example 1 

Consider as the first example the relatively simple optical im- 
age whose header is given in Table 4. The NAXIS and NAXIS j 
keywords indicate that we have a four-dimensional image con- 
sisting of 512 columns, 512 rows, 196 planes, and one polar- 
ization. The degenerate STOKES axis is used simply to convey 
a coordinate value which applies to every pixel in the image. 
The CRPIX j keywords tell us that the reference point is at pixel 
coordinate (256, 257, 1, 1). The VCiJa keywords default to the 
unit matrix which indicates that no bulk rotation or shear is ap- 
plied to the pixel coordinates. Intermediate world coordinates 
may thus be computed from Eq. (1) as 







( -0.003 








0) 




'Pi 


- 256 \ 


y 







0.003 










P2 


- 257 


z 










7128.3 







Pi 


- 1 


, s , 




, 








I, 




,P4 


- 1 J 



Since 'VELOCITY' and 'STOKES' are linear axis types, the 
velocity and Stokes value of each point are found simply by 



Table 5. Sample calculations for example 1. 



parameter 


units 


SE corner 


NE corner 


NW corner 




pixel 


(1,2) 


(1,512) 


(511,512) 


(P3,P4) 


pixel 


(1, 1) 


(1,1) 


(196, 1) 


X 


deg 


0°765000 


a 765000 


-0°765000 


y 


deg 


-0°765000 


0°765000 


0°765000 


</> 


deg 


45° 000000 


135° 000000 


225° 000000 


e 


deg 


88°918255 


88°918255 


88°9 18255 


a 


deg 


47°503264 


47°595581 


44°064419 


S 


deg 


62° 7951 11 


64° 324332 


64°324332 


Velocity 


ms~' 


500000.00 


500000.00 


1890018.50 


Stokes 




1.0 = 1 


1.0 = 1 


1.0 = 1 



adding the coordinate value at the reference point to the relative 
coordinate. Thus, 

Velocity = 500000.0 + 7128.3 fe-l)ms"', 
Stokes = 1 (I polarization) . 

The CTYPE1 and CTYPE2 keywords denote a TAN 
(gnomonic) projection for which (x,y) are projection plane co- 
ordinates for the zenithal perspective projection with the source 
of the projection at the center of the sphere. Thus the native lon- 
gitude and latitude are given by 

4> = arg(-y,x), 
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Fig. 35. Rgipro jection) - _/? fl (SIN) in arcsec plotted versus 9 for various projections. 



, ( 180° 1 1 
9 = tan" 1 ; , 

which, on substitution, become 

(j> = arg(/? 2 -257,/?i -256)+ 180°, 



lVCPi-256) 2 + Cp 2 -257) 2 J' 

where we have used Eqs. (15), (14), and (55). 

The celestial coordinate system is equatorial since the 
CTYPE ia begin with RA and DEC and the RADESYSa and 
EQUINOXa cards denote that these are in the IAU 1984 sys- 
tem. Zenithal projections have (<f>o, 9o) = (0, 90°) for which the 
CRVAL/ give equatorial coordinates a p = 45°83 in right ascen- 
sion and 6 P = 63°57 in declination. The equatorial north pole 
has a longitude of 180° in the native coordinate system from the 
LONPOLEa keyword. LATPOLEa is never required for zenithal 
projections and was not given. Thus, Eqs. (6) for the right as- 
cension and declination become 

sin 6 = sin 9 sin(63°57) - cos 9 cos cos(63°57) , 
cos 5 sin(or - 45°83) = cos 9 sin <p , 
cos£cos(a-45°83) = sin0cos(63°57) 

+ cos 9 cos (p sin(63°57) . 

Sample calculations for points on the diagonal near the three 
corners of the image are shown in Table 5. 

If we define the projection non-linearity as the departure of 
9 from r = ^x 1 +y 2 , then in this 1°5 x 1°5 image it amounts 
to ~ 075 at the corners. However, in a 6° x 6° image it quickly 
escalates to ~ 0.'5, sixty times larger. Comparison of a for the 
southeast and northeast corners indicates the significant effect 
of grid convergence even in this moderate sized image. The two 
differ by 22! 16, most of which is attributable to the cos 6 term 
in converting longitude offsets to true angular distances. The 
effect of grid convergence is small for (ao, So) near the equator 
and large near the poles. 



Figure 35 investigates the effect of projective non- 
linearities for moderate field sizes for the commonly used 
zenithal projections. It shows the difference in Re between the 
various projections and the SIN projection as a function of na- 
tive latitude 9. The SIN projection is used as reference since it 
always has R g less than the other zenithal projections. The dif- 
ference for all projections exceeds 1 arcsec for values of 9 less 
than 88° and the difference for the TAN projection exceeds one 
milliarcsec only 440 arcsec from the native pole. Such nonlin- 
earities would be significant in optical, VLBI, and other high 
resolution observations. 



7.3.2. Header interpretation example 2 

While the previous header was a realistic example it overlooked 
many of the concepts introduced in this paper. Consider now 
the header given in Table 6. Although contrived to illustrate as 
much as possible in one example, this header could conceivably 
arise as a "tile" from a conic equal area projection of a region 
of the southern galactic hemisphere centered at galactic coor- 
dinates {{, b) = (90°, -25°). The tile size of 2048 x 2048 pixels 
is approximately 10° x 10°, and this particular tile is situated 
immediately to the north of the central tile. 

The header defines ecliptic coordinates as an alternate co- 
ordinate system, perhaps to help define the distribution of zodi- 
acal light. This has the same reference point and transformation 
matrix as the primary description and the reference values are 
translated according to the prescription given in Sect. 2. The 
RADESYSA card indicates that the ecliptic coordinates are mean 
coordinates in the IAU 1984 system, though the author of the 
header has been sloppy in omitting the EQUINOXA card which 
therefore defaults to J2000.0. Note that, in accord with Paper I, 
all keywords for the alternate description are reproduced, even 
those which do not differ from the primary description. 
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Table 6. Second example FITS header (blank lines have been inserted for clarity). 



123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789 



NAXIS 






2 


/ 


2 -dimensional image 


NAXIS1 = 






2848 


/ 


x axis (fastest) 


NAXIS2 = 






2848 


/ 


y axis (2nd fastest) 


HJD-OBS = 




44258, 


,7845612 


/ 


MJD at start of observation 


CRPIX1 = 






1024.5 


/ 


Pixel coordinate of reference point 


CRPIX2 = 






-1823.5 


/ 


Pixel coordinate of reference point 


PC1_1 






1 


/ 


Transformation matrix element 


PC1_2 






-8.884 


/ 


Transformation matrix element 


PC2_1 






-8.882 


/ 


Transformation matrix element 


PC2_2 






1 




Transformation matrix element 


CDELT1 = 






-8.885 


/ 


x-scale 


CDELT2 = 






8.885 


/ 


y-scale 


CTYPE1 = 


'GLON- 


-COE' 




/ 


Conic equal area projection 


CTYPE2 = 


'GLAT- 


-COE' 




/ 


Conic equal area projection 


PV2_1 






-25.8 


/ 


Conic mid-latitude 


CRVAL1 = 






98.8 


/ 


Galactic longitude at reference point 


CRVAL2 = 






-25.8 


/ 


Galactic latitude at reference point 


CRPIX1A = 






1824.5 


/ 


Pixel coordinate of reference point 


CRPIX2A = 






-1823.5 


/ 


Pixel coordinate of reference point 


PC1_1A = 






1 


/ 


Transformation matrix element 


PC1_2A = 






-8.884 


/ 


Transformation matrix element 


PC2_1A = 






-8.882 


/ 


Transformation matrix element 


PC2_2A = 






1 


/ 


Transformation matrix element 


CDELT1A = 






-8.885 


/ 


x-scale 


CDELT2A = 






8.885 


/ 


y-scale 


TTVPF1 A — 








/ 
/ 


V-UIL1C trljUal died pi UJ CL L1UIL 


CTYPE2A = 


'ELAT- 


■COE' 




/ 


Conic equal area projection 


PV2_1A = 






-25.8 


/ 


Conic mid-latitude 


CRVAL1A = 




-7, 


,0308934 


/ 


Ecliptic longitude at reference point 


CRVAL2A = 




34. 


,8474143 


/ 


Ecliptic latitude at reference point 


LONPOLEA= 




6 


,3839706 


/ 


Native longitude of ecliptic pole 


LATPOLEA= 




29, 


,8114400 


/ 


Ecliptic latitude of native pole 


RADESYSA= 


'FK5 






/ 


Mean IAU 1984 ecliptic coordinates 



The problem will be to determine the galactic and eclip- 
tic coordinates corresponding to a field point with pixel coor- 
dinates (1957.2,775.4). We begin as before by substituting in 
Eq.(l) 




-0.005 0.00002 W/n - 1024.5 \ 
-0.00001 0.005 )\P2 + 1023.5 J ' 



Note that this transformation matrix describes a slight skewness 
and rotation; the x-, and >-axes correspond to the following two 
non-orthogonal lines in the pixel plane; 

p 2 = 0.002 (pi - 1024.5)- 1023.5, 
p 2 = 250 (pi - 1024.5) - 1023.5 . 

The values of the (x,y) projection plane coordinates for the 
chosen point are shown in Table 7 together with the remain- 
ing calculations for this example. 

The next step is to compute the native longitude and lati- 
tude. Like all standard conies, the conic equal area projection 
has two parameters, # a and rj, given by the keywords ?Vi_la 
and PV;'_2a attached to latitude coordinate ;'. In this example 



Table 7. Calculations for example 2. 



(PuPi) 


1957.2 


775.4 


(x,y) 


-4°6275220 


8°9851730 




-25°0000000 


0! 0000000 


61,82 


-25°0000000 


-25°0000000 


y 


-0.8452365 




C,¥ 


-0.4226183 


-122° 8711957 


(<P,ff) 


-4°7560186 


-15°8973800 




-90°0000000 


90! 0000000 


(f,b) 


85°2439814 


-15°8973800 


(^p,ySp) 


-179°9767827 


29°81 14400 




-14°7066741 


43° 0457292 



PV2_1 is given but the author has again been sloppy in omitting 
PV2_2, which thus defaults to 0. Explicit inclusion of this key- 
word in the header would obviate any suspicion that it had been 
accidentally omitted. The native coordinates, (0, 6) are obtained 
from (x,y) using Eqs. (117) and (129) via the intermediaries C, 
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Yq, and R g given by Eqs. (125), (127), and (116). These in turn 
are expressed in terms of y, 6\, and 82 given by Eqs. (128), 
(119), and (120). The calculations are shown in Table 7. 

At this stage we have deprojected the field point. In this 
example the alternate coordinate description defines the same 
projection as the primary description. Indeed, it may seem odd 
that the formalism even admits the possibility that they may 
differ. However, this is a realistic possibility, for example in 
defining multiple optical plate solutions based on the TAN pro- 
jection. It now remains to transform the native spherical coordi- 
nates into galactic coordinates, {(, b), and ecliptic coordinates, 
(A.,/3). To do this we need to apply Eqs. (2) and in order to do 
that we need (( p , b p ) and (A p ,fi p ). 

Looking first at galactic coordinates, we have ((o,bo) = 
(90°, -25°) and (p p = 0. This conic projection has (0o,#o) = 
(0°, -25°), and because bo = 8q and (p p = 0°, the native and 
galactic coordinate systems must coincide to within an off- 
set in longitude. However, it is not obvious what l p should 
be to produce this offset. Equation (8) has one valid solution, 
namely b v = 90°. The special case in point 2 in the usage 
notes for Eqs. (9) and (10) must therefore be used to obtain 
(£ p , b p ) = (-90°, 90°). The galactic coordinates of the field 
point listed in Table 7 are then readily obtained by application 
of Eqs. (2). 

The header says that the ecliptic coordinates of the refer- 
ence point are (A ,/3 Q ) = (-7°0300934, 34°8474143) and the 
native longitude of the ecliptic pole is <f> p = 6°3839706. It also 
specifies LATPOLEA as 29°8 1 14400. In this case Eq. (8) has two 
valid solutions, /3 P = -25! 1367794 ± 54^9482194, and the one 
closest in value to LATPOLEA (in fact equal to it) is chosen. If 
LATPOLEA had been omitted from the header its default value 
of +90° would have selected the northerly solution anyway, 
but of course it is good practice to make the choice clear. The 
value of A p may be obtained by a straightforward application of 
Eqs. (9) and (10), and the ecliptic coordinates of the field point 
computed via Eqs. (2) are listed in Table 7 as the final step of 
the calculation. The reader may verify the calculation by trans- 
forming the computed galactic coordinates of the field point to 
mean ecliptic coordinates. 

7.3.3. Binary table representations of example 2 

Table 8 shows the FITS header for a set of images stored in 
binary table image array column format. The images are similar 
to that described in header interpretation example 2, Sect. 7.3.2, 
with primary image header illustrated in Table 6. 

In this example the images are stored as 2-D arrays in 
column 5 of the table and each row of the table contains a 
2048 x 2048 pixel image of a different region on the sky. This 
might represent a set of smaller images extracted from a sin- 
gle larger image. In this case all coordinate system parameters 
except for the reference pixel coordinate are the same for each 
image and are given as header keywords. The reference pixel 
coordinate for the primary and secondary description are given 
in columns 1 to 4 of the binary table. 

Table 9 shows the header for the same image given in pixel 
list format. There are 10000 rows in this table, each one listing 



the pixel coordinates (XPOS, YPOS) of every detected "event" 
or photon in the image. For illustration purposes, this table 
also contains an optional 'DATA_QUALITY' column that could 
be used to flag the reliability or statistical significance of each 
event. A real image may be constructed from this virtual image 
as the 2-dimensional histogram of the number of events that 
occur within each pixel. The additional TLMINn and TLMAXn 
keywords shown here are used to specify the minimum and 
maximum legal values in XPOS and YPOS columns and are use- 
ful for determining the range of each axis when constructing 
the image histogram. 

7.3.4. Header interpretation example 3 

This example has been adapted from a real-life FITS data file. 
The simplicity of the header shown in Table 10 is deceptive; it 
is actually presented as an example of how not to write a FITS 
header, although the latent problem with its interpretation is 
quite subtle. 

Observe that the image spans 180° in native longitude and 
90° in native latitude and that the reference pixel lies outside 
the image. In fact, the reference pixel is located so that the na- 
tive longitude runs from 45° to 225° and hence the image lies 
partly inside and partly outside the normal range of native lon- 
gitude, [-180°, 180°]. 

In fact, as might be expected, this makes no difference 
to the computation of celestial coordinates. For example, in 
computing the celestial coordinates of pixel (1,1) we readily 
find from Eqs. (83) and (84) that the native coordinates are 
(0, 0) = (225°, -45°). The fact that 4> exceeds 180° becomes ir- 
relevant once Eqs. (2) are applied since the trigonometric func- 
tions do not distinguish between <p = 225° and (p = -135°. 
The latter value is the appropriate one to use if the cylinder of 
projection is considered to be "rolled out" over multiple cycles. 
Consequently the correct galactic coordinates are obtained. 

The problem only surfaces when we come to draw a coordi- 
nate grid on the image. A meridian of longitude, for example, is 
traced by computing the pixel coordinates for each of a succes- 
sion of latitudes along the segment of the meridian that crosses 
the image. As usual, in computing pixel coordinates, the ce- 
lestial coordinates are first converted to native coordinates by 
applying Eqs. (5), and the native longitude will be returned in 
the normal range [-180°, 180°]. Consequently, in those parts of 
the image where (f> > 180° the pixel coordinates computed 
will correspond to the point at - 360°, i.e. in the part of the 
principle cycle of the cylindrical projection outside the image. 

In principle it is possible to account for this, at least in spe- 
cific cases, particularly for the cylindrical projections which 
are somewhat unusual in this regard. In practice, however, it 
is difficult to devise a general solution, especially when similar 
problems may arise for projection types where it is not desir- 
able to track (p outside its normal range. For example, consider 
the case where a Hammer- Aitoff projection is used to represent 
the whole sky; since its boundary is curved there will be out- 
of-bounds areas in the corner of the image. Normally a grid 
drawing routine can detect these by checking whether the in- 
verse projection equations return a value for $ outside its nor- 



34 Mark R. Calabretta and Eric W. Greisen: Representations of celestial coordinates in FITS 

Table 8. Example binary table image array header (blank lines have been inserted for clarity). 



123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789 
XTENSION= 'BINTABLE' / Binary table extension 



BITPIX 






8 


/ 


8-bit bytes 


NAXIS 


_ 




2 


/ 


2-dimensional binary table 


NAXIS1 


_ 


16777232 


/ 


Width of table in bytes 


NAXIS2 


_ 




4 


/ 


Number of rows in table 


PCOUNT 


_ 







/ 


Size of special data area 


GCOUNT 


_ 




1 


/ 


One data group (required keyword) 


TFIELDS 


= 




5 


/ 


number of fields in each row 


TTYPE1 


= 


' 1CRP5 ' 




/ 


Axis 1: reference pixel coordinate 


TF0RM1 


= 


'1J 




/ 


Data format of column 1: 1*4 integer 


TTYPE2 


= 


' 2CRP5 




/ 


Axis 2: reference pixel coordinate 


TF0RM2 




'1J 




/ 


Data format of column 2 : 1*4 integer 


TTYPE3 




' 1CRP5A ' 




/ 


Axis 1A: reference pixel coordinate 


TF0RM3 




■U 




/ 


Data format of column 3 : 1*4 integer 


TTYPE4 




' 2CRP5A ' 




/ 


Axis 2A: reference pixel coordinate 


TF0RM4 




'1J 




/ 


Data format of column 4: 1*4 integer 


TTYPE5 




' Image ' 




/ 


2-D image array 


TF0RM5 




'4194304J' 




/ 


Data format of column 5 : 1*4 vector 


TDIM5 




'(2048,2048) 




/ 


Dimension of column 5 array 


MJD0B5 




44258 


.7845612 


/ 


MJD at start of observation 


COMMENT 


The following keywords define the primary coordinate description 


COMMENT 


of the images 


contained 


in Column 5 of the table. 


11PC5 






1 


/ 


Transformation matrix element 


12PC5 






-0.004 


/ 


Transformation matrix element 


21PC5 






-0.002 


/ 


Transformation matrix element 


22PC5 






1 


/ 


Transformation matrix element 


1CDE5 






-0.005 


/ 


Axis 1: scale 


2CDE5 






0.005 


/ 


Axis 2: scale 


1CTY5 




'GLON-COE' 




/ 


Axis 1: conic equal area projection 


2CTY5 




'GLAT-COE' 




/ 


Axis 2: conic equal area projection 


2PV5_1 






-25.0 


/ 


Conic mid-latitude 


1CRV5 






90.0 


/ 


Axis 1: galactic longitude at reference point 


2CRV5 






-25.0 


/ 


Axis 2: galactic latitude at reference point 


COMMENT 


The following keywords define the secondary coordinate description 


COMMENT 


of the images 


contained 


in Column 5 of the table. 


11PC5A 






1 


/ 


Transformation matrix element 


12PC5A 






-0.004 


/ 


Transformation matrix element 


21PC5A 






-0.002 


/ 


Transformation matrix element 


22PC5A 






1 


/ 


Transformation matrix element 


1CDE5A 






-0.005 


/ 


Axis 1A: scale 


2CDE5A 






0.005 


/ 


Axis 2A: scale 


1CTY5A 




'ELON-COE' 




/ 


Axis 1A: conic equal area projection 


2CTY5A 




'ELAT-COE' 




/ 


Axis 2A: conic equal area projection 


2PV5_1A 






-25.0 


/ 


Conic mid-latitude 


1CRV5A 




-7 


.0300934 


/ 


Axis 1A: ecliptic longitude at reference point 


2CRV5A 




34 


.8474143 


/ 


Axis 2A: ecliptic latitude at reference point 


L0NP5A 




6 


.3839706 


/ 


Native longitude of ecliptic pole 


LATP5A 




29 


.8114400 


/ 


Ecliptic latitude of native pole 


RADE5A 




'FK5 




/ 


Mean IAU 1984 ecliptic coordinates 


EQUI5A 






2000.0 


/ 


Coordinate epoch 



END 
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Table 9. Example pixel list header (blank lines have been inserted for clarity). 



123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789 



XTENSION= 


' BINTABLE ' 




/ 


Binary table extension 


BITPIX 


= 




8 


/ 


8-bit bytes 


NAXIS 


= 




2 


/ 


2-dimensional binary table 


NAXIS1 


= 




5 


/ 


Width of table in bytes 


NAXIS2 


= 




10000 


/ 


Number of rows in table 


PCOUNT 


= 







/ 


Size of special data area 


GCOUNT 


= 




1 


/ 


One data group (required keyword) 


TFIELDS 


= 




3 


/ 


Number of fields in each row 


TTYPE1 


= 


' DATA_QUALITY ' 


/ 


Quality flag value of the photon 


TF0RM1 


= 


'IB 




/ 


Data format of the field: 1-byte integer 


TTYPE2 




'XPOS 




/ 


Axis 1: pixel coordinate of the photon 


TF0RM2 




'11 




/ 


Data format of column 2 : 1*2 integer 


"T"T MTITI 

TLHIN2 






1 


/ 


Lower limit of axis in column 2 


TLMAX2 






2048 


/ 


Upper limit of axis in column 2 


TTYPE3 




'YPOS 




/ 


Axis 2: pixel coordinate of the photon 


TF0RM3 




'11 




/ 


Data format of column 3: 1*2 integer 


TLMIN3 






1 


/ 


Lower limit of axis in column 3 


TLMAX3 






2848 


/ 


Upper limit of axis in column 3 


MJD0B3 




44258, 


,7845612 


/ 


MJD at start of observation 


COMMENT 


The following keywords define the primary coordinate description. 


TCRP2 






1024.5 


/ 


Axis 1: reference pixel coordinate 


TCRP3 






-1823.5 


/ 


Axis 2: reference point pixel coordinate 


TPC2_2 






1 


/ 


Transformation matrix element 


TPC2_3 






-8.884 


/ 


Transformation matrix element 


TPC3_2 






-8.882 


/ 


Transformation matrix element 


TPC3_3 






1 


/ 


Transformation matrix element 


TCDE2 






-8.885 


/ 


Axis 1: scale 


TCDE3 






8.885 


/ 


Axis 2: scale 


TCTY2 




'GLON-COE' 




/ 


Axis 1: conic equal area projection 


TCTY3 




'GLAT-COE' 




/ 


Axis 2: conic equal area projection 


TPV3_1 






-25.8 


/ 


Conic mid-latitude 


TCRV2 






98.8 


/ 


Axis 1: galactic longitude at reference point 


TCRV3 






-25.8 


/ 


Axis 2: galactic latitude at reference point 


COMMENT 


The following keywords define the secondary coordinate description. 


TCRP2A 






1824.5 


/ 


Axis 1A: reference pixel coordinate 


TCRP3A 






-1823.5 


/ 


Axis 2A: reference point pixel coordinate 


TP2_2A 






1 


/ 


Transformation matrix element 


TP2_3A 






-8.884 


/ 


Transformation matrix element 


TP3_2A 






-8.882 


/ 


Transformation matrix element 


TP3_3A 






1 


/ 


Transformation matrix element 


TCDE2A 






-8.885 


/ 


Axis 1A: scale 


TCDE3A 






8.885 


/ 


Axis 2A: scale 


TCTY2A 




'ELON-COE' 




/ 


Axis 1A: conic equal area projection 


TCTY3A 




'ELAT-COE' 




/ 


Axis 2A: conic equal area projection 


TV3_1A 






-25.8 


/ 


Conic mid-latitude 


TCRV2A 




-7, 


,0308934 


/ 


Axis 1A: ecliptic longitude at reference point 


TCRV3A 




34, 


,8474143 


/ 


Axis 2A: ecliptic latitude at reference point 


L0NP3A 




6. 


,3839706 


/ 


Native longitude of ecliptic pole 


LATP3A 




29, 


,8114400 


/ 


Ecliptic latitude of native pole 


RADE3A 




' FK5 




/ 


Mean IAU 1984 ecliptic coordinates 


EQUI3A 






2000.0 


/ 


Coordinate epoch 



END 
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Table 10. Third example FITS header. 



123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789 



NAXIS 


= 


2 / 2 -dimensional image 


NAXIS1 


= 


181 / x axis (fastest) 


NAXIS2 


= 


91 / y axis (2nd fastest) 






*? f\ fil / Pt vol /" , /"\f*T"'/*l"ir"i^i"t"*3 r\f t 1 a f qt 1 a~n f* a Y\r\ \ t"i +■ 
L* D . w / rlACl LUUiUlIldLc UI Iclcl CILLc pu±IlL 


CRPIX2 




46.8 / Pixel coordinate of reference point 


CDELT1 




-1.8 / x-scale 


CDELT2 




1.8 / y-scale 


CTYPE1 


= 'GLON-CAR' 


/ Plate carree projection 


CTYPE2 


= 'GLAT-CAR' 


/ Plate carree projection 


CRVAL1 




38.8 / Galactic longitude at reference point 


CRVAL2 




35.8 / Galactic latitude at reference point 



mal range. It may thus determine the proper boundary of the 
projection and deal with the discontinuity that arises when a 
grid line passes through it. 

How then should the header have been written? Note that 
the problem exists at the lowest level of the coordinate descrip- 
tion, in the conversion between (x,y) and (<t>,6), and the so- 
lution must be found at this level. The problem arose from a 
particular property of cylindrical projections in that they have 
x oc (/>. We must use this same property, which we might call 
"^-translation similarity", to recast the coordinate description 
into a more manageable form, ^-translation similarity simply 
means that changing the origin of (f> corresponds to shifting the 
image in the x-direction. In other words, we can transfer the 
reference point of the projection from its current location to 
another location along the native equator without having to re- 
grid the image. The fact that the VCiJa matrix is unity in this 
example makes this task a little simpler than otherwise. 

Note first that because the image straddles (p = 180° we 
can't simply reset CRPIX1 so as to shift the reference point to 
4> = -360°; the image would then straddle (p = -180°, which 
is no improvement. In this example it is convenient and suffi- 
cient to shift the reference point to {<p,9) = (180°, 0°), which 
corresponds to pixel coordinate {p\,pi) = (46.0,46.0). Hence 
we need to reset CRPIX1 to 46.0 and adjust the keywords which 
define the celestial coordinate system. The reader may readily 
verify that the galactic coordinates of the new reference point 
are {I, b) = (210°, -35°) and whereas the old, implied value of 
LONPOLE was 0° when So > 0, now that So < its new im- 
plied value is 180°, and this is correct. However, we will set it 
explicitly anyway. The keywords to be changed are therefore 

CRPIX1 = 46.0, 
CRVAL1 = 210°, 
CRVAL2 = -35°, 
LONPOLE = 180°. 

What if the PC i.ja matrix was not unity? The problem of 
determining the pixel coordinates where ((j), 9) = (-180°, 0°) 
would have presented little extra difficulty, although in general 
CRPIX2 would also need to be changed. On reflection it may 
come as a surprise that changing CRPIX j like this does not 
fundamentally alter the linear transformation. However it may 



readily be verified that the only effect of changing CRPIX j is to 
change the origin of the (x,y) coordinates. 

7.4. Header construction examples 

This section considers the more difficult problem faced by FITS 
writers; that of constructing world coordinate system headers. 

Example 1 is contrived to illustrate the general meth- 
ods used in constructing a spherical coordinate representation. 
Paying homage to Claudius Ptolemy, it actually constructs a 
terrestrial coordinate grid for the Mediterranean region as seen 
from space. 

Example 2 constructs headers for the infra-red dust maps 
produced by Schlegel et al. (1998) who regridded data from 
the COBE/DIRBE and IRAS/ISSA surveys onto two zenithal 
equal area projections. 

Example 3 considers the case of long-slit optical spec- 
troscopy. It is concerned in particular with solving the problem 
of producing three world coordinate elements for a data array 
of only two dimensions. 

Example 4 constructs a dual coordinate description for an 
image of the moon and considers the problem of producing 
consistent scales for each. It also suggests an extension to deal 
with the rings of the planet Saturn. 

7.4.1. Header construction example 1 

Our first example of FITS header construction concerns satel- 
lite photography of the Earth. Figure 36 shows a part of 
the world which probably would have been recognizable to 
Claudius Ptolemy. 

A satellite 2230 km immediately above Cairo aims its dig- 
ital camera directly towards Athens and adjusts the orientation 
and focal length to include Cairo in the field near the bottom 
edge of the frame. The 2048 x 2048 pixel CCD detector ar- 
ray employed by the camera is centered on its optical axis so 
Athens is at pixel coordinate (1024.5, 1024.5). Cairo - at the 
satellite's nadir - is later found to be at (681.67,60.12). The 
task is to overlay a terrestrial coordinate grid on this image. 

For the sake of simplicity we will assume pin-hole cam- 
era optics. Figure 37 identifies the geometry as that of a tilted, 
near-sided zenithal perspective projection with the nadir at the 
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Fig. 36. The Earth from 2230 km above Cairo looking directly towards 
Athens. 



Fig. 37. The geometry of Fig. 36. Cairo at C is at the reference point 
of the projection with Athens at A. The camera with aperture at P and 
focal plane F is not drawn to scale, though the rest of the diagram is. 
Note that ji < so -fi, as marked, is a positive value. 



reference point. The point of projection P corresponds to the 
pin-hole of the camera and the tilted plane of projection is par- 
allel to the camera's focal plane F. The image in the focal plane 
is simply scaled (and inverted) with respect to that on the plane 
of projection. Clearly a tilted AZP projection is a better match 
to the geometry than SZP in this instance. 

The satellite altitude of 2230km is 0.350 Earth radii and 
since the projection is near-sided we may immediately write 
^ = -1.350. 

Determination of y requires knowledge of the coordinates 
of Cairo (31°15E,30°03N) and Athens (23°44E,38°00N). 
From spherical trigonometry we may deduce that the angu- 
lar separation between the two cities is 10°2072 and also that 
Athens is on a bearing 36°6252 west of north from Cairo. 

Now since Cairo is at the native pole of the projection the 
angular separation between the two cities is just their difference 
in native latitude, 90° - 9a- Using the sine rule in triangle PAO 
in Fig. 37 we obtain 

90 o -6» = -y-sin- 1 0usinr). (194) 

This equation, which takes account of the fact that fi is negative, 
may be solved iteratively to obtain y = 25° 845 8. 

The obliqueness of this view of the Earth, occasioned by 
the fact that the image is not oriented along a meridian, is han- 
dled partly via LONPOLEa and partly as a bulk image rotation 
via PC/_/'a. Note that these are not complementary; they pro- 
duce distinct effects since the tilted AZP projection does not 
have point symmetry. Figure 37 shows the situation - the gen- 
erating sphere with Cairo at the native pole must be oriented so 
that Athens is at native longitude 180°. The native longitude of 
the terrestrial pole (substituting for the celestial pole in this ex- 
ample) is offset from this by the bearing of Athens from Cairo, 
i.e. 4> p = 180° - 36°6252 = 143°3748. 

Since Athens is at native longitude 180° its x-coordinate, 
like Cairo's, must be zero. The inequality of the corresponding 



pixel coordinates must have arisen by the satellite rotating the 
camera about its optical axis thereby rotating the CCD detec- 
tor in the focal plane. The angle of this bulk rotation is read- 
ily deduced from the pixel coordinates of Athens and Cairo, 
arg(j A - jcJk - z'c) = 19°570. This rotation is to be applied 
via PCLja. The direction of rotation is completely determined 
by the requirement that Athens' x-coordinate be zero, regard- 
less of the sign of CDELT/a, and the resulting PCLja rotation 
matrix is shown below. 

The reference pixel coordinates, CRPIX ja, were mea- 
sured from the image and the reference coordinates for Cairo, 
CRVAL/a, were obtained from an atlas, so the last remaining 
unknowns are the scales, CDELT/a. From previous calculations 
we know that Athens is at (<f>, 0) = (180°, 79°7928) and we may 
apply Eqs. (20) and (21) to obtain (x,y) = (0,8°7424). The 
distance in pixel coordinates between the two cities is readily 
found to be 1023.5 so the y-scale must be 8°7424/1023.5 = 
0°008542 per pixel. 

The x-scale cannot be determined like this since both cities 
have the same x-coordinate. However, the x-, and y-scales must 
be equal because the focal plane is parallel to the plane of pro- 
jection and ray-tracing through the pin-hole therefore results in 
an isotropic change of scale. Do not confuse this with the state- 
ment made in Sect. 5.1.1 that with y + the projection is not 
scaled true at the reference point; this refers to the differential 
scale between (x,y) and (</>, 6). 

Though the image in the focal plane is inverted through 
the pin-hole, thus indicating a negative scale, we can assume 
that the camera compensated by reading out the CCD array in 
reverse order. Thus for the image of Fig. 36 we make both of 
the CDELT/a positive in order to have east to the right as befits a 
sphere seen from the outside. Putting this all together we have 

NAXIS = 2, 
NAXIS1 = 2048, 
NAXIS 2 = 2048, 
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CRPIX1 




681 .67 , 


CRPIX2 


= 


60.12, 


PC1_1 


= 


0.9422 , 


PC1_2 


_ 


-0.3350, 


PC2_1 


= 


0.3350, 


PC2_2 


= 


0.9422 , 


CDELT1 


_ 


0°008542 , 


CDELT2 


_ 


0°008542 , 


CTYPE1 


= 


'TLON-AZP' , 


CTYPE2 


= 


'TLAT-AZP' , 


PV2_1 




-1.350, 


PV2_2 




25°8458, 


CRVAL 1 




31°15 , 


CRVAL2 




30°03 , 


LONPOLE 




143°3748 , 


WCSNAME 




'Terrestrial coordinates' 



This example was simplified by the fact that the nadir was 
known and was included in the field of view. This probably 
would not be the case in a more realistic example where there 
may also be some uncertainty in the target coordinates and the 
value of p. In such cases the mapping parameters would have 
to be determined by least squares via the identification of an 
adequate number of surface features with known coordinates. 

7.4.2. Header construction example 2 

The following example comes from the 4096 x 4096 pixel in- 
frared dust maps produced by Schlegel, Finkbeiner, & Davis 
(SFD, 1998). The authors chose to regrid data from the 
COBE/DIRBE and IRAS/ISSA maps onto two zenithal equal 
area (ZEA) projections centered on the galactic poles. The pro- 
jection formula given in their Appendix C expressed in terms 
of standard 1 -relative FITS pixel coordinates (pupi) is 

pi - 1 = 2048 Vl -nsinb cos I + 2047.5 , 
p 2 - 1 = -«2048 Vl - nsinb sin t + 2047.5 , 

where n — +1 for the north Galactic pole and n = -1 for the 
south Galactic pole maps. Now for ZEA from Eqs. (1), (12), 
(13), and (69) 

180° , 

CDELTl(/?i - CRPIX1) = V2 Vl-sin0sin0, 

7T 

180° 

CDELT2 (p2 - CRPIX2) = - V2 Vl -sin6»cos0. 

n 

If we take the north Galactic pole case first, the SFD equations 
may be rewritten as 

pi - 2048.5 = -2048 Vl -sinfe sin(£ - 90°) , 
p 2 - 2048.5 = -2048 Vl - sin b cos({ - 90°) . 

By inspection of the two sets of equations we must have 

NAXIS = 2, 



NAXIS1 


= 4096, 


NAXIS2 


= 4096, 


CRPIX1 


= 2048.5, 


CRPIX2 


= 2048.5, 


CDELT1 


= -180° V2/ (2048?!-) , 


CDELT2 


= 180° V2/ (20487T), 


CTYPE1 


= 'GLON-ZEA' , 


CTYPE2 


= 'GLAT-ZEA' , 


I 


= + 90°, 


b 


= e. 



Now, writing I and b in place of a and 6 in Eqs. (2), we have 
to determine f p , b p , and <f> p to give ({, b) in terms of (</>, 6). This 
is easy because we know b p = 90° and the simple special case 
Eqs. (3) apply, so 

{ = 4> + ({ p ~4> p - 180°). 

We have a degree of freedom here since only ( p - <p p is deter- 
mined. It's best to let (p p take its default value of 0° for zenithal 
projections with the celestial pole at the native pole (it's 180° 
otherwise), so we must have 

CRVAL 1 = 270°, 
CRVAL 2 = 90°, 
LONPOLE = 0°. 

Although LONPOLE assumes its default value here it would of 
course be wise to write it explicitly into the header. The proce- 
dure for the south Galactic pole case is similar. The SFD equa- 
tions may be rewritten 

Pi - 2048.5 = 2048 Vl - sin(-fo) sin(90° - €) , 
p 2 - 2048.5 = 2048 Vl - sin(-fe) cos(90° - €) , 
whence 

CRPIX1 = 2048.5, 
CRPIX2 = 2048.5, 
CDELT1 = 180° V2/(2048tt), 
CDELT2 = -180° V2/(2048tt), 

t = 90° - 4> , 

b = -9. 

Equations (4) apply for b p = -90° so 

e = (^p + P )-0. 

Again we let (p p take its default value of 180°, so 

CRVAL 1 = 270°, 
CRVAL 2 = -90°, 
LONPOLE = 180°. 

It's generally easier to interpret coordinate headers than to 
construct them so it's essential after formulating the header to 
test it at a few points and make sure that it works as expected. 
Note that this sort of translation exercise wouldn't be necessary 
if the formalism of this paper was used right from the start, i.e. 
in the regridding operation used to produce the maps. 
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7.4.3. Header construction example 3 

Consider now the coordinate description for the two- 
dimensional image formed by a long slit spectrograph. We as- 
sume that the wavelength axis of length 1024 and dispersion 
AA nm/pixel corresponds to the p\ pixel coordinate, and the 
2048 pixel spatial axis corresponds to p2- The slit is centered 
on equatorial coordinates («o, <5o) and oriented at position angle 
p measured such that when p = the first spectrum is north- 
wards. We will assume that the telescope and spectrograph op- 
tics are such that the distance along the slit is in proportion to 
true angular distance on the sky with a separation between pix- 
els of cr arcsec. We do not consider curvature of the slit here, 
that is a distortion of the sort to be handled by the methods of 
Paper IV. 

Equiscaling along the length of the slit together with the 
one-dimensional nature of the spatial geometry indicate that 
any projection could be used that is equiscaled along a great 
circle projected as a straight line. Many projections satisfy this 
criterion. However, in practice the zenithal equidistant (ARC) 
projection is the most convenient choice. 

The one-dimensional spatial geometry may at first seem 
problematical, with each point along the slit having a different 
(a, 6). However, this is easily handled by introducing a third, 
degenerate axis and introducing a rotation. The rotation in posi- 
tion angle may be handled via the linear transformation matrix. 
Moreover, since we've opted for a zenithal projection, bearing 
in mind the discussion of Sect. 7. 1 , the rotation can also be han- 
dled via (p p (i.e. L0NP0LE). We will demonstrate both methods. 

Use of p is perhaps more straightforward, and the header 
may be written without further ado: 



NAXIS 




3, 


NAXIS1 




1024, 


NAXIS2 




2048, 


NAXIS 3 




1, 


CRPIX1 




1, 


CRPIX2 




1024.5 , 


CRPIX3 




1, 


CDELT1 




AA, 


CDELT2 




-cr/3600, 


CDELT3 




1, 


CTYPE1 




'WAVELEN ' , 


CTYPE2 




'RA ARC 


CTYPE3 




'DEC- -ARC 


CUNIT1 




'nm' , 


CRVAL1 




A , 


CRVAL2 




«o , 


CRVAL3 




*0, 


LONPOLE 




90° +p. 



LONPOLE is the only card which requires explanation. Since 
no rotation is introduced by the linear transformation, the slit, 
which coincides with the /?2-axis, maps directly to the x-axis. 
However, because CDELT2 is negative, the two axes run in op- 



posite directions and, given that when p = the /?2-axis runs 
from north to south, x must run from south to north. Referring 
to the left side of Fig. 3 for zenithal projections, we see that 
when p = 0, the north celestial pole must be at = 90°. Since 
position angle increases in an anticlockwise direction (i.e. north 
through east) in the plane of the sky, the celestial pole must ro- 
tate clockwise as p increases. Thus, as p increases so must <p p , 
hence p = 90° + p. 

To verify this representation we will test it with (ao, <5o) = 
(150°,-35°),p = 30°, andcr = 2" for pixel coordinate ( 1 , 1, 1). 
The corresponding (s, x,y) coordinates are (0, 0°56861 1 1, 0°). 
The wavelength is therefore Aq. From Eqs. (15), (14), and (67) 
the native longitude and latitude are (<p, 6) = (90°, 89°43 13889), 
is always ±90°, and the value of 6 corresponds to an off- 
set of 2047" from the center of the slit as it should. From 
Eqs. (2) the right ascension and declination are (a, 5) = 
(150°3450039, -34°5070794), which are in the correct quad- 
rant. 

Most optical telescopes are better described by the TAN pro- 
jection. In this case the only difference in the header would be 

CTYPE2 = 'RA TAN' , 

CTYPE3 = 'DEC- -TAN' . 

Verifying it with the same values as before yields the same 
(s, x,y) and thus the same wavelength. From Eqs. (15), (14), 
and (54) we get (<f>,6) = (90°, 89°43 14076), the offset of 
2046'.'933 differing by only 67 mas. The right ascension and 
declination are (a, 6) = (150°3449926, -34°5070956). 

As an alternative, the original header could also have been 
written with the CTYPE i interchanged, so that the slit, still coin- 
cident with the ;?2-axis, maps directly to the y-axis. The header 
would be as above but with the following changes: 

CTYPE2 = 'DEC- -ARC , 

CTYPE 3 = ' RA ARC ' , 

CRVAL2 = (J , 
CRVAL3 = a , 
LONPOLE = 180° +p. 

LONPOLE differs because the slit now runs along the y-axis 
rather than the x-axis; the negative sign on CDELT2 is still 
needed to make y run counter to p2- 

There are several practical ways to rewrite the header using 
the PCLj or CDLj matrices. Looking first at the CDLj form we 
have 





' CDl 


.1 


CD1_2 


CD1_3 1 




[Pi 


- CRPIX1 ' 


x - 


= CD 2 


.1 


CD2_2 


CD2_3 




P2 


- CRPIX2 


,y , 


v CD 3 


.1 


CD3_2 


CD3_3 t 




,P3 


- CRPIX3 , 



Since the third axis is degenerate, with CRPIX3 = 1, we have 
p 3 - CRPIX3 = 0, so the values of CD1_3, CD2_3, and CD3_3 
are irrelevant. Moreover, CD1_2, CD2_1, and CD3_1 are all zero, 
hence 

s = CD1_1 (pi - CRPIX1), 

x = CD2_2 (p 2 - CRPIX2) , (195) 
y = CD3_2 (p 2 - CRPIX2) . 
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Effectively this provides (x, y) coordinates along the slit via the 
parametric equations of a line, where pi - CRPIX2 is the dis- 
tance parameter. 

As before, the main problem is to determine the correct an- 
gle of rotation. For zenithal projections <f> p = 180° by default, 
and referring to the left side of Fig. 3 we see that this corre- 
sponds to the y-axis. Thus when p = we require a rotation of 
90° to transform the /?2-axis onto the y-axis. For p > we need 
to rotate further so the angle of rotation is 90° + p. 

Therefore, CDELT; and LONPOLE in the original header 
must be replaced with 

CD1_1 = AA, 

CD2_2 = cos(90° + p)(<x/3600) , 
CD3_2 = - sin(90° + p)(<x/3600) , 
CD2_3 = 1, 
CD3_3 = 1. 

The negative sign on CD3_2 is associated with the rotation, not 
the scale. Note that, although the value of CD3_3 is irrelevant, 
it must be set non-zero otherwise the zero defaults for CDLj 
would make the third column zero, thereby producing a singu- 
lar matrix. Likewise, in this example and similarly below, we 
set CD2_3 non-zero to prevent the second row of the matrix be- 
coming zero for values of p such that cos(90° + p) = 0. 

The PC/_/ matrix formulation is similar but allows more 
flexibility in the way the scale is handled: 

1. Since CDELT / defaults to unity, omitting it allows PCL/ to 
duplicate the functionality of CD/_j. However, since PC/_/ 
defaults to the unit matrix rather than zero, there is no need 
to set PC3_3 specifically, hence 

PC 1.1 = AA, 

PC2_2 = cos(90° + p)(cr/3600) , 
PC3_2 = - sin(90° + p)(<x/3600) , 
PC2_3 = 1 . 

2. Equation (195), as a simple scaling relation, suggests the 
use of CDELT/. However, PC3_2 must be non-zero since x 
and y both depend on p2- Setting it to unity and allowing 
the remaining PC Lj to default to the unit matrix we have 

PC3_2 = 1, 

CDELT 1 = AA, 

CDELT2 = cos(90° + p)(<x/3600) , 

CDELT3 = - sin(90° + p)(<x/3600) . 

3. The "orthodox" method is to encode the rotation and scale 
separately in PC/_/ and CDELT/: 

PC2_2 = cos(90°+p), 

PC3_2 = -sin(90°+p), 

PC2_3 = 1, 

CDELT 1 = AA, 

CDELT2 = cr/3600, 

CDELT 3 = cr/3600. 



Arguably this is more amenable to human interpretation, 
especially if thoughtful comments are added. 

The above six methods should all be regarded as equally le- 
gitimate. In fact, there are infinitely many ways to encode this 
header, though most would disguise the essential simplicity of 
the geometry. 

7.4.4. Header construction example 4 

Thompson (1999) has applied the methods of this paper to the 
definition of solar coordinates for a variety of coordinate sys- 
tems. As the final header construction example we will con- 
sider a specific coordinate description for another solar system 
body, the Moon. 

A short exposure plate taken at Sydney Observatory at 
1:00 am on 15 February 1957 AEST (GMT +1000) of the full 
Moon near lunar perigee has been digitized and converted to a 
4096 x 4096 pixel image in FITS format. The scale is 1" per 
pixel with the moon centered in the image, and it is desired 
to construct a dual coordinate description. The first system is 
geocentric apparent equatorial coordinates in a gnomonic pro- 
jection. The second is selenographic coordinates in a zenithal 
perspective projection attached to the surface of the Moon. 

The ephemeris gives the geocentric apparent right ascen- 
sion of the Moon at the time as 09 h 41 m 13 ! il and declination as 
+08°34'26". Diurnal parallax would have caused the Moon to 
appear slightly offset from this position as seen from Sydney 
Observatory, but to a good approximation the coordinate sys- 
tems may be defined with the center of the Moon at the stated 
geocentric coordinates. The image was digitized in the normal 
orientation with north up and east to the left so the header for 
the primary description is straightforward: 



NAXIS 




2, 


NAXIS1 




4096, 


NAXIS2 




4096, 


MJD-OBS 




35883.625 , 


CRPIX1 




2048.5 , 


CRPIX2 




2048.5 , 


CDELT 1 




-0.0002778 , 


CDELT2 




0.0002778 , 


CTYPE1 




'RA TAN' 


CTYPE2 




'DEC- -TAN' 


CRVAL1 




145.30458 , 


CRVAL2 




8.57386, 


LONPOLE 




180.0, 


RADESYS 




' GAPPT ' . 



If any rotation had been introduced it could have been corrected 
for in the linear transformation matrix. 

The first step in constructing the secondary description is 
to determine the distance of the observer from the Moon. The 
ephemeris gives the equatorial horizontal parallax as 61'2973, 
corresponding to a distance between the centers of the Earth 
and Moon of 55.91 Earth radii. However, the observer may be 
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Fig. 38. Geometry of header construction example 4 (not to scale). 
The observer is at O, a distance of Moon radii from the center of 
the Moon. The Moon is projected as a near-sided zenithal perspective 
projection (pi < - 1), and the celestial sphere as a gnomonic projection. 
An alternative plane of projection is shown. 



closer or further away by up to one Earth radius depending on 
location and time of day and this 2% effect is deemed worthy of 
correction. At Sydney Observatory (longitude 10 h 04 m 49 s .2, lat- 
itude -33°51'41") the Greenwich apparent sidereal time was 
00 h 51 m 43 s .6. Thus the apparent right ascension and declina- 
tion of the zenith was 10 h 56 m 32 s :8, -33°51'41". The vector dot 
product then gives the distance correction as 0.69 Earth radii. 
Using the ratio of the Earth and Moon radii of 3.670 the cor- 
rected distance of 55.22 Earth radii indicates that the distance 
parameter for the zenithal perspective projection is fi = 202.64 
Moon radii. 

We now need to determine the correct relative scale. 
Figure 38 shows the geometry of the two projections where we 
note, by analogy with Fig. 4, that fi < -1. From the diagram 
we have 

^tan-U^). (196) 

The figure shows that this equation is not influenced by the 
choice of the plane of projection. Evaluating the derivative 
d/31 Ay at y = we find the relative scale is l/(|ju| - 1). 

The ephemeris records that the selenographic coordinates 
of the Earth at the time were {(, b) = (+0°26, +6°A5) and the 
position angle of the Moon's axis was C = 19°03. However, 
since the Earth subtends an angle of over 2° in the lunar sky, 
topocentric optical libration - the correction for the observa- 
tory's location - is significant. Application of the correction 
formulae derived by Atkinson (1951) gives the selenographic 
coordinates of Sydney Observatory as ((', V) = (-0°23, 5!89) 
and C = 18°94. Since the image was taken in the normal 



orientation and we have a zenithal projection it is convenient 
to account for the position angle by setting p = 180° - C. 
Adopting keyword values of SELN and SELT for selenographic 
coordinates we may write 



CRPIX1S 




2048.5 , 


CRPIX2S 


= 


2048.5 , 


CDELT1S 


= 


0.0562896 , 


CDELT2S 




0.0562896 , 


CTYPE1S 




'SELN-AZP' , 


CTYPE2S 




'SELT-AZP' , 


PV2_1S 




202.64 , 


CRVAL1S 




-0.23 , 


CRVAL2S 




5.89, 


LONPOLES 




161.06, 


WCSNAMES 




' SELENOGRAPHIC COORDINATES ' 



We have used the letter 'S' to denote the alternate coordinate 
system simply to demonstrate that there is no requirement to 
start with ' A' . A WC SN AME a keyword, denned inPaperI,isused 
to identify the coordinate system. Note that selenographic co- 
ordinates are right-handed so that selenographic longitude in- 
creases towards the west on the celestial sphere as seen from 
the Earth and this is handled by setting CDELT1S positive. 

As an extension of this example, a FITS header with 
three coordinate systems might be constructed for an image of 
Saturn; a celestial grid for the background, a saturnographic 
system for the surface of the planet, and a third system for its 
rings. The rings might be described by a zenithal equidistant 
(ARC) projection with associated linear transformation matrix 
set to match the oblique viewing angle. 

7.5. Realization 

Calabretta (1995) has written, and made available under a GNU 
license, a package of routines, wcslib, which implements all 
projections and coordinate conversions defined here. It contains 
independent C and Fortran libraries. 

The Fortran library includes a routine, pgsbox, which is 
based on pgplot and uses wcslib to draw general curvilinear 
coordinate axes. It was used to generate Fig. 36. 

8. Summary 

We have developed a flexible method for associating celestial 
coordinates with a FITS image and implemented it for all pro- 
jections likely to be of use in astronomy. It should be a rel- 
atively simple matter to add new projections should the need 
ever arise. 

The FITS-header keywords defined in this paper are sum- 
marized in Table 1 1 with the 3 -letter projection codes listed in 
Table 12. The column labeled 6 Q gives the native latitude of the 
reference point in degrees (<f>o = always) and the parameters 
associated with each projection are listed in the nomenclature 
of Sect. 5. 
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Table 11. Summary of new coordinate keywords introduced in this paper; see also Table 2 for alternate types used in binary tables. 



Keyword Type Sect. Use 



Status 



Comments 



LONPOLEa floating 2.2 coordinate rotation new 



LATPOLEa floating 2.4 coordinate rotation new 



RADESYSa character 3.1 



frame of reference new 



EQUINOX a floating 3.1 coordinate epoch new 



EPOCH floating 3.1 

MJD-OBS floating 3.1 



coordinate epoch 
time of observation 



deprecated 
new 



Longitude in the native coordinate system of the 

celestial system's north pole. 

Default: 0° if S > 6 , 180° otherwise. 

Latitude in the native coordinate system of the 

celestial system's north pole, or equivalently, the 

latitude in the celestial coordinate system of the 

native system's north pole. 

Default: 90°, unless (0„, 6 , (p P - <f>o) = (0, 0, ±90°) 

in which case it has no default. 

Reference frame of equatorial and ecliptic 

coordinates; recognized values are given in Table 1. 

Default: 'FK4' if EQUINOXa < 1984.0, 'FK5' if 

> 1984.0, or ICRS if EQUINOXa is not given. 

Epoch of the mean equator and equinox in years; 

Besselian if RADESYSa is FK4 or FK4-N0-E, 

Julian if FK5; not applicable to ICRS or GAPPT. 

Default: EPOCH if given, else 1950.0 if RADESYSa is FK4 

or FK4-N0-E, or 2000.0 if FK5. 

Replaced by EQUINOXa. 

Modified Julian Date (JD - 2400000.5) of observation. 
Default: DATE-OBS if given, else no default. 
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Appendix A: Mathematical methods 

A.1. Coordinate rotation with matrices 

The coordinate rotations represented in Eqs. (2) or (5) may be 
represented by a matrix multiplication of a vector of direction 
cosines. The matrix and its inverse (which is simply the trans- 
pose) may be precomputed and applied repetitively to a variety 
of coordinates, improving performance. Thus, we have 



'M (r n 
m = r 2 \ 
, n ) [r 3 i 


rn rn\( V " 
r 22 r 23 m' , 
r 32 r 33 ){n' / 


where 




(l',m',n) = 


(cos S cos a, cos 5 sin a, sin 5) , 


(l,m,ri) 


(cos 6 cos </>, cos sin (/>, sin 6) , 


• = 


- sin o-p sin <p p - cos a p cos (f> p sin 6 P 


m = 


cos a p sin p - sin a p cos <f> p sin 6 p , 


rn = 


cos 4> p cos 6 P , 


r 2 \ = 


sin ap cos (p p - cos a p sin (p p sin 6 p , 


r 22 = 


- cos a p cos (p p - sin a p sin (p p sin 6 P 


rn = 


sin p cos 6 P , 


rn = 


cos a p cos (S p , 


r 32 = 


sin a p cos 6 P , 


r 33 = 


sin <5 P . 



The inverse equation is 







' rn r 2 \ r 31 N 




' I ^ 


m' 




rn r 22 r 32 




m 


, n' , 




, rn r 23 r 33 t 




, n , 
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FITS 

Code 


00 


00 


Projection name 


Projection parameters associated with latitude 1 axis i 
PVL8a PVLla PVL2a PVL3a PVLma 


A7P 


u 


QO° 


Zenithal perspective 


H y 


C7D 


u 


Qfl° 


Slant zenithal perspective 


ii /h Q 


TAN 


0° 


90° 


Gnomonic 




STG 


0° 


90° 


Stereographic 




SIN 


0° 


90° 


Slant orthographic 


£ n 


ARC 


0° 


90° 


Zenithal equidistant 




ZPN 


0° 


90° 


Zenithal polynomial 


Po Pi Pi P3 P m form = 0,. 


ZEA 


0° 


90° 


Zenithal equal-area 




AIR 


0° 


90° 


Airy 


ft. 


CYP 


0° 


0° 


Cylindrical perspective 




CEA 


0° 


0° 


Cylindrical equal area 


A 


CAR 


0° 


0° 


Plate carree 




MER 


0° 


0° 


Mercator 




SFL 


0° 


0° 


Sanson-Flamsteed 




PAR 


0° 


0° 


Parabolic 




MOL 


0° 


0° 


Mollweide 




AIT 


0° 


0° 


Hfitntnpr- A i tnfF 

1 llll 1 11 1 1L 1 ill LUll 




COP 


0° 


#a 


Conic perspective 




COE 


0° 




f^nnir pmifll-arpa 




COD 


0° 


On 


Conic equidistant 


u d 'I 


COO 


0° 


a 


Conic orthomorphic 


ft, n 


BON 


0° 


0° 


Bonne's equal area 




PCO 


0° 


0° 


Polyconic 




TSC 


0° 


0° 


Tangential Spherical Cube 




CSC 


0° 


0° 


COBE Quadrilateralized Spherical Cube 




QSC 


0° 


0° 


Quadrilateralized Spherical Cube 





f Parameters PVLOa, PVLla, and PVL2a associated with longitude axis implement user-specified values of (<f>o,0o) as discussed in Sect. 2.5, 
and PVL3a, and PVL4a encapsulate the values of LONPOLEa and LATPOLEa respectively, as discussed in Sect. 2.6. 



A. 2. Iterative solution 

Iterative methods are required for the inversion of several of 
the projections described in this paper. One, Mollweide's, even 
requires solution of a transcendental equation for the forward 
equations. However, these do not give rise to any particular dif- 
ficulties. 

On the other hand, it sometimes happens that one pixel and 
one celestial coordinate element is known and it is required to 
find the others; this typically arises when plotting graticules on 
image displays. Although analytical solutions exist for a few 
special cases, iterative methods must be used in the general 
case. If, say, P2 and a are known, one would compute pixel co- 
ordinate as a function of 6 and determine the value which gave 
P2- The unknown pixel coordinate elements would be obtained 
in the process. 

This prescription glosses over many complications, how- 
ever. All bounded projections may give rise to discontinuities 
in the graph of pi versus S (to continue the above exam- 
ple), for example where the meridian through a crosses the 
(f> = ±180° boundary in cylindrical, conic and other projec- 
tions. Even worse, if the meridian traverses a pole represented 
as a finite line segment then p2 may become multivalued at a 
particular value of 6. The derivative dp2/dS will also usually 



be discontinuous at the point of discontinuity, and it should 
be remembered that some projections such as the quad-cubes 
may also have discontinuous derivatives at points within their 
boundaries. 

We will not attempt to resolve these difficulties here but 
simply note that wcslib (Calabretta, 1995) implements a solu- 
tion. 

Appendix B: The slant orthographic projection 

The slant orthographic or generalized SIN projection derives 
from the basic interferometer equation (e.g. Thompson et al., 
1986). The phase term in the Fourier exponent is 

p = (e-eo)-B, (B.l) 

where eo and e are unit vectors pointing towards the field center 
and a point in the field, B is the baseline vector, and we measure 
the phase p in rotations so that we don't need to carry factors 
of 2n. We can write 

p = p u u + p v v + p w w , (B.2) 

where (u, v, w) are components of the baseline vector in a right- 
handed coordinate system with the w-axis pointing from the 
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Table A.l. Projection aliases. 



Name 



Alias for 



Name 



Alias for 



Gnomonic 

= Central 
Near-sided perspective 
Clarke's (first) 
Clarke's (second) 
James' 
La Hire's 

Approximate equidistant 

zenithal 
Approximate equal area 

zenithal 
Postel 

= Equidistant 

= Globular 
Lambert azimuthal 
equivalent 

= Lambert azimuthal 
equal area 

= Lambert polar 
azimuthal 

= Lorgna 
Gall's cylindrical 
Cylindrical equal area 
Simple cylindrical 

= Central cylindrical 

= Cylindrical central 
perspective 

= Gall's stereographic 
Lambert's cylindrical 

= Lambert's equal area 
Behrmann equal area 
Gall's orthographic 

= Approximate Peter's 
Lambert's equal area 



AZP with 11 = 

AZP withyu = 1.35 
AZP withyu = 1.35 
AZP withyu = 1.65 
AZP withyu = 1.367 
AZP withyu = 1.71 

AZP withyu = 1.7519 

AZP with n = 2 AIM 
ARC 



ZEA 



CYP withyu = l,A = V2/2 
CYP with fi = 00 
CYP with n = 0, A = 1 



CYP with// = oo,i=l 

CEA with A = 3/4 
CEA with,! = 1/2 

CEA with A = 1 



Miller 

Equirectangular 
Cartesian 

= Equidistant 
cylindrical 

Cassini 

Transverse Mercator 
= Transverse cylindrical 
orthomorphic 

Sinusoidal 

= Global sinusoid (GLS) 
= Mercator equal-area 
= Mercator-Sanson 
= Sanson's 

Craster 

Bartholomew's atlantis 
Mollweide's homolographic 

= Homolographic 

= Homalographic 

= Babinet 

= Elliptical 
Hammer equal area 

= Aitoff 

= Aitov 
Bartholomew's nordic 
One-standard conic 

= Tangent conic 
Two-standard conic 

= Secant conic 
Murdoch conic 
Alber's 

= Alber's equal area 
Lambert equal area 
Lambert conformal conic 
Werner's 



CAR with x scaled by 2/n 
CAR with unequal scaling 
CAR 



CAR transverse case 
MER transverse case 



SFL 



PAR 

MOL oblique case 
MOL 



AIT 



AIT oblique case 
Conic with 8\ = 6*2 

Conic with ft * 6 2 

similar to COD 
COE 

COE with 2 = 90° 

for spherical Earth = COO 

BON with ft = 90° 



geocenter towards the source and the M-axis lying in the equa- 
torial plane, and 

p u = cos 6 sin <p , 

p v = - cosflcos^, (B.3) 
p w = sin 6 - 1 , 

are the coordinates of (e - eo), where (<p, 0) are the longitude 
and latitude of e in the spherical coordinate system with the 
pole towards eo and origin of longitude towards negative v, as 
required by Fig. 3. Now, for a planar array we may write 

n u u + n v v + n w w = (B.4) 

where (n u , n v , n w ) are the direction cosines of the normal to the 
plane. Using this to eliminate w from Eq. (B.2) we have 

P = [Pu~ —Pw\u + [p v - —p w ]v (B.5) 

flw 1l\v 



Being the Fourier conjugate variables, the quantities in brackets 
become the Cartesian coordinates, in radians, in the plane of the 
synthesized map. Writing 

f = njn w , 
J] = n v /n w , 

Eqs. (61) and (62) are then readily derived from Eqs. (B.3) and 
(B.5). For 12-hour synthesis by an east-west interferometer the 
baselines all lie in the Earth's equatorial plane whence , rf) = 
(0, cottfo), where So is the declination of the field center. For a 
"snapshot" observation by an array such as the VLA, Cornwell 
& Perley (1992) give (£,77) = (-tanZsin^, tanZcos^f), where 
Z is the zenith angle and x is the parallactic angle of the field 
center at the time of the observation. 

In synthesizing a map a phase shift may be applied to the 
visibility data in order to translate the field center. If the shift 
applied is 

Ap = q u u + q v v + q w w (B.7) 
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where (q u , q v , q w ) is constant then equation (B.2) becomes 

p = (p u - q u )u + (p v - q v )v + (p w - q w )w , (B.8) 

whence equation (B.5) becomes 



(B.9) 



[(Pv - qv) - n(p w - q w )]v. 

Equations (61) and (62) become 

180° r , 180° r 

x = [ cos 6 sin <f> + £ (sin 6 - 1) J [q u - gq w \ , 

n n 

180° r ,. lx , 180° r 

y = [ cos #cos (p -rj (sin 8 - 1) J [q v - rjq w \ , 

n ' n 

(B.10) 

from which on comparison with Eqs. (61) and (62) we see that 
the field center is shifted by 

180° 

(Ax, Ay) = (q u - gq w , q v - rjq w ) . (B.l 1) 

7T 

The shift is applied to the coordinate reference pixel. 



Appendix C: Projection aliases 

Table A. 1 provides a list of aliases which have been used by 
cartographers for special cases of the projections described in 
Sect. 5. 
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